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Abstract 


Combined heat and power plants (CHPs) enable an efficient co-generation of heat and 
electricity and thus are crucial for a resource efficient energy supply. As the share of 
renewable generation is increasing in electric grids, the operation of CHPs gets more 
and more challenging. The traditional heat-driven operation of CHPs is not suitable 
to react to the volatile availability of renewable energy, to variable electric demand as 
well as to fluctuating electricity prices. However, flexible operation of CHPs requires 
thermal storage. As an alternative to dedicated thermal storage tanks, the water in- 
side a heating grid can be used as thermal storage if grid dynamics are considered 
in operations planning. It is very complex to find good or, ideally, optimal opera- 
tion schedules for CHPs when the dynamics of a heating grid are included. These 
dynamics are dominated by a variable-dependent time delay, which results in a non- 
convex problem formulation. Hence, common optimization approaches reach their 
limits and do not find the global optimum when the thermal dynamics of a heating 
grid are considered. This thesis presents the following new methods to neverthe- 
less find optimal schedules for operation of CHPs considering heating grid dynamics: 
i) the first method solving problems with variable-dependent time delays to global 
optimality by proposing a novel outer approximation of the pipe outflow tempera- 
ture in a primal-dual global optimization scheme, ii) a method introducing a hybrid 
discrete-continuous time grid and discretized temperatures to enable an accurate rep- 
resentation of the variable-dependent time delays in a mixed-integer linear program 
and iii) an approach allowing a measurement based identification of the heating grid 
dynamics that enables an easy application to real world grids. 


Kurzfassung 


Kraft-Wärme-Kopplung (KWK) ermöglicht eine effiziente Erzeugung von Wärme und 
Elektrizität und spielt somit eine zentrale Rolle für eine ressourcenschonende En- 
ergieversorgung. Mit stetig wachsender Einspeisung erneuerbarer Energieträger ins 
Stromnetz steigt die Volatilität am Strommarkt und regelbare Erzeuger müssen flex- 
ibel auf diese Schwankungen reagieren können. Die traditionell wärmegeführte Be- 
triebsweise von KWK-Anlagen ist deshalb nicht mehr zeitgemäß. Für einen flexibleren 
Betrieb von KWK-Anlagen sind jedoch Wärmespeicher notwendig. Nicht nur in Spe- 
ichertanks, sondern auch im Wärmenetz selbst kann Wärmeenergie gespeichert wer- 
den. Hierzu muss die Dynamik des Wärmenetzes in der Einsatzplanung der KWK- 
Anlagen berücksichtigt werden. Die Modellierung und Optimierung dieser thermis- 
chen Dynamik ist jedoch komplex, da sie auf einer variablen-abhängigen Verzögerung 
beruht. Diese variablen-abhängige Verzögerung führt zu einer nicht-konvexen Prob- 
lemstellung, für die herkömmliche Optimierungsverfahren nicht ausgelegt sind. Um 
trotzdem optimale Lösungen zu erreichen, werden in dieser Arbeit folgende neue 
Methoden zur optimalen Einsatzplanung von KWK-Anlagen unter Berücksichtigung 
der Dynamik von Wärmenetzen entwickelt: i) die erste Methode, die Probleme mit 
variablen-abhängigen Verzögerungen global optimal löst, indem eine neue äußere Ab- 
schätzung der Rohraustrittstemperatur als duales Problem in einem iterativen, glob- 
alen Optimierungalgorithmus eingesetzt wird, ii) eine Methode, die mit Hilfe einer 
Diskretisierung der möglichen Vorlauftemperaturen sowie einem hybriden diskret- 
kontinuierlichen Zeitstrahl eine exakte Formulierung der Verzögerung als gemischt- 
ganzzahliges Problem erlaubt und iii) eine Methode, die eine messwertgestützte Mod- 
ellidentifikation der Wärmenetzdynamik ermöglicht und somit einfach in realen Net- 
zen eingesetzt werden kann. 
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FD volume at temperature level / consumed in time slot t and leaving pipe 
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Symbol Description 
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k value of digit 

P upper bound of position of digit 
p lower bound of position of digit 
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1 Introduction 


The reduction of green house gas emissions and an efficient use of energy resources 
are a major concern in today’s societies. District heating grids play an important role 
for a sustainable energy supply as they enable efficient co-generation of heat and elec- 
tricity in combined heat and power plants (CHPs) as well as the use of thermal energy 
normally dissipated to the environment e.g. from industrial production processes or 
waste incineration. District heating grids transport thermal energy from heat sources 
to heat sinks using mostly water as transport medium. This water is heated at gen- 
eration sites to a supply temperature chosen by the heating grid operators. Then the 
hot water is transported via supply pipes to heat consumers. They cool down the 
hot water to a return temperature and the water is transported back to the generation 
sites via separate return pipes. 


Worldwide there are about 80,000 heating grid systems supplying about 14 EJ ther- 
mal energy per year of which 6,000 are located in Europe supplying about 2.5 EJ 
thermal energy each year [Wer17]. The operators of these heating grids are facing 
multiple challenges today. On the one hand, they need to supply the thermal demand 
as resource-, energy- and cost-efficiently as possible. On the other hand, energy and 
especially electricity markets are becoming more and more volatile due to a rising 
share of renewable power generation. Hence, the traditional heat-driven operation 
strategy for CHPs reaches its limits and operators need to consider intra-day fluctua- 
tions of prices and availability of electricity. Thus, a closely coupled coordination of 
co-generation of different energy forms is needed, which is one of the central topics 
of 4th generation district heating [LWW 14]. 


In particular, a better integration of different energy sectors provides benefits for re- 
newable power integration into electric grids. In electric power systems, generation 
needs to equal demand at all times. As renewable generation is not fully controllable, 
energy storage is of major importance to balance demand and supply in power grids 
with a high share of renewable generation. Storage units that are able to store en- 
ergy for minutes, hours, days and months are required to stabilize these future grids. 
Pumped-hydro storage, the most affordable form of electric energy storage, as well 
as thermal energy storage tanks are both very suitable for storage cycles of up to a 
few days. However, investment costs for thermal energy storage tanks are orders of 
magnitude cheaper than for pumped-hydro storage [L@C*16]. Hence, using ther- 
mal energy storage instead of electric storage reduces investment costs importantly. 
In addition, thermal energy can directly be stored in a heating grid for a few hours 


2 1 Introduction 


thanks to its dynamic behavior [Gro12]. Thus, every heating grid offers thermal stor- 
age capacity without any investment costs, if the thermal dynamics of the grid are 
considered in operations planning. 


However, using this inherent grid storage in operations planning is a very complex 
task, as the thermal dynamics of a heating grid are driven by a mass-flow dependent 
transport time, which leads to a time delay in the arrival of an increase of temperature 
from the beginning to the end of a pipe. When applying optimization approaches to 
operations planning problems for district heating grids, this variable-dependent delay 
results in a non-convex problem formulation. Due to its importance for future energy 
systems, several publications have addressed this operational optimization of heating 
grids. As shown in Chapter 2, iterative and sequential solution approaches have been 
proposed to cope with the non-convex variable-dependent time delay. However, none 
of these approaches prove global optimality of their solution. Thus, it is unknown if 
they leverage the full potential of the heating grid as thermal energy storage capacity. 
In addition, the algorithms that have been proposed in the literature are only rarely 
used in real world applications as they are complicated to implement and use. 


Therefore, this work addresses the following research questions: 


e Can we develop an algorithm which reliably finds the global optimal solution for 
operational optimization of combined heat and power plants whilst considering the 
dynamics of heating grids? 


e Can we find exact and easy to engineer algorithms which offer fast solution times 
and thus are applicable for real world heating grids? 


To tackle the first question, Chapter 3 presents “multiparametric delay modeling” the 
first approach reaching global optimality for operational planning of heating grids. 
“Multiparametric delay modeling” introduces binary variables that represent the cur- 
rent transport delay and are used to derive an outer-approximation of the outflow 
temperature of a pipe. Chapter 4 and Chapter 5 address the second question, facilitat- 
ing optimization of heating grids in real-world operations planning. The hybrid-time 
grid approach presented in Chapter 4 focuses on a fast to solve but still very accurate 
model. It adapts optimization models recently developed for pipeline scheduling in 
the petrochemical industry by introducing discrete temperature levels and a hybrid 
discrete-continuous time grid. The main focus of the delay matrix approach presented 
in Chapter 5 lies on an easy and fast parameterization of the optimization model. The 
model dynamics are represented in a static delay matrix which can be calculated 
using historic temperature and flow measurements. Chapter 6 then compares the op- 
timization approaches presented in this thesis to each other as well as with a literature 
approach using small example heating grids. In addition, a real-world case study is 
presented that demonstrates the scalability and benefits of the approach of Chapter 5. 
Finally, Chapter 7 summarizes and discusses the results of this work. 


2 State of the art 


As mentioned in the introduction, several past publications address the optimiza- 
tion of heating grids. This chapter gives first an overview on the state of the art for 
operational optimization of heating grids. Second, theory on global optimization of 
non-convex problems is introduced with focus on non-convex problem classes that 
are required to describe thermal dynamics of heating grids. In this literature review 
it becomes apparent that current optimization theory does not provide suitable meth- 
ods for global optimization of heating grid dynamics. Third, Section 2.3 introduces 
pipeline scheduling for crude oil products which models similar flow dynamics as 
found in heating grids. 


2.1 Operational optimization of district heating grids 


Trying to find optimal set points for operation of district heating grids has been stud- 
ied in many past publications. Common goals in operational optimization of district 
heating grids are the reduction of green house gas emissions as well as the reduc- 
tion of operational cost. However, different optimization problems aiming at these 
goals can be distinguished in literature: i) If multiple heat production units are part 
of the grid, a key question is to find their schedules determining which unit produces 
when [SP08]. ii) Finding the optimal supply temperature is of high interest, which 
ensures security of supply for all consumers while reducing thermal and hydraulic 
losses [AMP14, VTD17]. A combination of both problems is possible, resulting in a 
fully dynamic optimization of the heating grid [BBR95]. 


In this thesis, the main focus lies on scheduling of production units. However, the 
global optimization approach presented in Chapter 3 offers a fully dynamic optimiza- 
tion. 


2.1.1 Scheduling of CHPs without consideration 
of thermal dynamics 


For scheduling of CHPs, district heating grid dynamics are very often neglected and 
a static heat load forecast is assumed. With these assumptions, it is common to use 
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a mixed-integer linear programming (MILP)! discrete time unit commitment model 
like, e.g., in [Dot97, NNM00, ROGJ17]. This section introduces such a basic discrete 
time MILP model for scheduling of CHPs without consideration of thermal dynamics. 
Existing applications and extensions of this model are reviewed at the end of this 
section. 


With the discrete time representation we get a set S+ of time slots t of equal duration. 
As input parameters do not vary within each time slot, model accuracy is not reduced. 
In addition, variables are constant within every time slot t. All CHPs are collected in 
set Sọ. For every CHP i in each time slot t, we introduce non-negative real variables 
representing the electric power output p;; and the thermal power output q;,; as well as 
the running status of the unit n; being a binary variable (0: off / 1: on). Lower limits 
of power and heat output are PŁ and Q} and the respective upper limits of power and 
heat output are PU and QU. To enforce that the CHP only produces power and heat 
when it is running (n;; = 1) we get 


n PE < pip < nip PU Vi € Sg, Yt € St, (2.1) 
nit QE < qip < nip QU Vi € Sg, VE € Sp. (2.2) 


Electric power output p;; and heat power output q;+ of a CHP are not independent. 
This relationship is influenced by the used generation technology. Gas turbines as 
well as gas engines usually have a linear connection between heat and electric output 
and hence a fixed power-to-heat ratio, whereas for CHPs with extraction condensing 
turbine the power-to-heat ratio is flexible and limited by a capability region as shown 
in Figure 2.1 [Cer02]. The feasible region of a CHP i with extraction condensing 
turbine shown in Figure 2.1 can be described with inequalities 


Pit Saı — Bı Git VEE Si, 
Pit 2 %2 — Ba Git Vt E Si, (2.3) 
Pit 2 —a3 + P3 ` Git VEE Si 


with a1, a2, &3, b1, p2 and 63 being parameters describing slope and offset of the 
limiting lines. 


Assuming that electric demand pferd and thermal demand gfe@nd are known, en- 


ergy balances for heat and electricity for every time slot t can be formulated. Because 
CHPs are usually connected to a large electricity grid, it is assumed that electric en- 
ergy can be sold and purchased. The volume of electricity sold in time slot t is p°", 
whereas the volume of electricity bought in time slot t is pit , both being non-negative 
real variables. If no energy storage like battery or pumped-hydro storage is installed, 


1 MILPs are optimization problems with real and integer (discrete) variables. The problem formulation 


consists of a linear objective function to be minimized or maximized as well as linear equality and 
inequality constraints linking real and integer variables. In many cases the integer variables are binary 
taking only two values. There are many off-the-shelf solvers to solve MILPs such as [BEN05, Gur16, 
FRV* 18]. 
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Figure 2.1: Feasible operation region of a CHP i with extraction condensing turbine linking normalized 
electric output p;+/ pH and normalized thermal output q;,:/ Qu following [Cer02]. 


electric energy needs to be produced at the same time it is consumed. Thus, the 
electric power balance becomes 


päemand _ Z Pis + py = piel! YEE Sp. (2.4) 
i€Sg 


For thermal energy there are usually no open markets enabling hourly sales and 
purchases, leading to 


german + gi = i gis VEE Sp. (2.5) 
i€Sg 
In this thermal energy balance, heat losses q/°** of the heating grid are accounted for 
as well. 


In unbundled electricity markets (like all electricity markets in Europe) utilities pro- 
ducing electricity do not need to consider transmission losses in the electric grid, as 
this is the responsibility of grid operators. They purchase energy at the energy mar- 
kets to compensate losses in transmission and distribution grids. Thus, for operations 
planning of a power producer grid losses are not relevant and do not appear in the 
electric power balance in (2.4). 


To complete the MILP formulation for scheduling of CHPs, the objective function 
is the last missing piece. As mentioned in Section 2.1, major goals for operation of 
district heating grids are the reduction of CO, emissions and the reduction of cost. If 
CO; emissions are considered in the generation cost, both goals can be combined by 
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minimizing the objective function 


min \ (price $ py ie price’!!! . pi") 
tes} 
+} (ag Niet eu ` Pit + ee. ait) . (2.6) 
LES ES, 
Here price?" and price?!" are time varying prices for buying and selling electricity 
and Cf"s! is the cost for running CHP i in one time slot. C? and C!**'#@" are the 
variable cost for the production of electricity (el) and heat of CHP i. 


Existing applications and extensions of the basic CHP scheduling model 


Many publications use a scheduling model as introduced above. For example, in 
[Dot97] the author proposes an economic dispatch model for CHPs which allows 
flexible generation by extending the model with thermal storage tanks. Time vary- 
ing cost of generation and time varying prices for sales to the energy markets are 
considered in scheduling of decentralized CHPs in [NNMO0]. As well considering 
decentralized CHPs with thermal storage tanks, experiences from a real-world on- 
line optimization are shared in [FKBV13]. Here, an intra-day optimization reacts to 
unforeseen variations in generation or load. Another work investigates how future 
electricity prices can influence the operations strategy of CHPs [ROGJ17]. A review 
of different methods for short term operations planning of CHPs without considera- 
tion of grid dynamics is given in [SP08]. Being not intended for large heating grids, 
the algorithm proposed in [LM17] as well does not consider time delays in tempera- 
ture propagation. It uses a Taylor series expansion around a given working point to 
linearize the remaining non-linear heating grid dynamics. 


2.1.2 Thermal dynamics of heating grids: the node method 


The model and the approaches presented in the previous section do not consider grid 
dynamics. To consider dynamics of a heating grid in simulation or optimization, 
the node method [BBR95] is a very common approach of modeling. Besides many 
academic works (e.g. [SFT*05, LWSt 16, ZFZC18]), commercial simulation tools rely 
on it [S0106] and simulations show that it achieves a better performance than other 
modelling approaches [Tr699]. 


Thus, this section introduces the basic principles of thermal dynamics of heating grids 
using the node method according to [BBR95]. As this and other work [Gro12], this 
thesis relies on the following assumptions: 


e incompressible transport medium 
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no leakages 

e perfect mixing in nodes (one single outflow temperature per node) 

e heat propagates only with the mass flow (advection), no conduction, no diffusion 
¢ return temperatures are constant 

e pressure dynamics are not considered as they are too fast 


e specific heat capacity and density of water are constant in the relevant temperature 
range (50 °C to 130 °C) 


e thermal demand of consumers is known and given for each time slot 


The key concept of the node method is to track the propagation of a volume of water 
at a certain temperature in the pipe. A district heating grid is described by nodes 
and pipes connecting the nodes. Thermal energy can be injected or withdrawn from 
the grid at nodes, if thermal generation sites or thermal loads are present. Besides 
generator parameters and load profiles, additional technical information like heat ca- 
pacities and pipe diameters is required. The node method is a discrete time approach 
calculating temperatures at and mass flows between the nodes. Based on past tem- 
peratures of nodes and past mass flows between the nodes, the node temperatures for 
the next time step are calculated. Temperatures at intermediate points between the 
nodes are not calculated. 


With the aforementioned assumptions (incompressible transport medium and no leak- 
ages) and assuming no storage tanks in nodes, there is a balance of mass flows m; 
flowing into and out of a node: 


L Mit = L Mit Vt = St. (2.7) 


icin flow icout flow 
Similarly, the mass flow into and out of a pipe are equal. 


Assuming perfect mixing of inflow temperatures T” results in a single outflow tem- 
perature Be for a node, which can be calculated as 


Licin flow (7 i itis) 


Lieinflow Mit 


ToS Vt © S. (2.8) 


Thermal generation injecting heat into the system consumes water from a node of 
the return pipe network, heats it and feeds it to a node of the supply pipe network. 
Thermal loads withdrawing heat from the system consume water from a node in the 
supply pipe network, chill it and feed it to a node of the return pipe network. Heating 
and chilling of water leads to a change of its inner energy Q; which can be described 
by 

Qi = Cpr Ci z TE VEE Sı (2.9) 


8 2 State of the art 


where ri; is the mass flow between supply and return side and c, is the specific heat 


capacity of the transport medium water (assumed to be constant). T; "" Y and men 
are the temperatures at supply and return side. 


loss,suppl 
pump . t Qe 


m; 
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= ) eneration load El 
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Figure 2.2: A heating grid with one generation and one load. Red lines: supply pipes; blue lines: return 
pipes; red arrows: thermal energy injection, retrieval and losses. 


Figure 2.2 shows the most simple setup for a heating network with one thermal gen- 
eration asset and one thermal load. Thermal losses occur in both supply and return 
network. However, thermal losses in the supply network are more important than 
in the return network due to the larger temperature difference to the surrounding 
area. Thermal loads usually have a controller adjusting the mass flow from supply 
to return side such that a given return temperature is reached for the arriving supply 
temperature and the current demand. As there is a pump close to the generation site 
with a constant pressure output, a valve at the load is sufficient for this mass flow 
control. The supply temperature is controlled at thermal generation sites. The return 
temperature depends on the return temperature arriving from the consumers. For 
all but one generation unit, the mass flow is controlled to fit the scheduled thermal 
energy supply of this generation unit. For the remaining generation unit, the mass 
flow and heat output remain variable to cope with unforeseen changes in thermal 
demand. 


With this very common control strategy, the operators of a heating grid can control the 
temperature levels only in the supply lines directly. The temperature level in return 
pipes mainly depends on the controllers at the consumers whose behavior is hard to 
predict. The ability to control the temperature is crucial for storing thermal energy 
in a heating grid [Gro12]. Therefore, only the thermal dynamics in the supply lines, 
which are under direct control of the heating grid operator, are considered in this 
thesis and a constant return temperature is assumed. 


As pressure dynamics travel with speed of sound in water (about 1200 m/s), they are 
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a lot faster than the temperature dynamics which travel with flow velocity (at most 
1.5 m/s to 2 m/s). Hence, pressure dynamics are neglected and only temperature 
dynamics are explained in the following. 


In the node method the outflow temperature T? “ut ofa pipe is calculated with a linear 
combination of past and present inflow temperatures T}". Example 2.1 shows how 
this linear combination represents the propagation of temperature in a pipe. 


Example 2.1: 


Figure 2.3 shows the temperature distribution inside a pipe of length L at the beginning of 
a time slot t. 


Figure 2.3: Scheme of pipe filled with different temperatures at beginning of a time slot t. The red- 
dashed volume leaves the pipe in time slot t. The blue arrow indicates the flow direction 
(left to right). 


T”, T?” , and 27 are the inflow temperatures feeding the pipe in past time slots t — 1, 
t—2and t — 3. Axı, Ax;_1, Axı_n and Ax;_3 are the distance traveled by the water inside 
the pipe in time slots t, t — 1, t — 2 and t — 3, respectively. 

As in time slot t the water moves by Ax, the red-dashed volume in Figure 2.3 leaves the pipe 
in time slot t. Hence, the average outflow temperature of this time slot T?"' is the average 
temperature of the red-dashed volume. If heat losses are not considered, the example shown 
in Figure 2.3 leads to (2.10) mixing the two temperatures inside the red-dashed volume 
with a weighted average: 


Di a Tas 
ree = 311-2 + z Ti-a (2.10) 


The propagation of the transport medium in the pipe Ax; in time slot t can be cal- 
culated using velocity v; or mass flow mų with At being the duration of a time slot, 
p being the density of the transport medium (water) and d being the diameter of the 
pipe. The resulting equation is 
At , 
Ax; = At -v = LLL Vt © St. (2.11) 
d 

er (3) 
Following Example 2.1, the pipe outflow temperature er in every time slot t can 
be calculated as a linear combination of past inflow temperatures T;” with maximum 
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max 


transport delay t’”* and weight factors a: 


max 


T . 
T Yo wets YEE Sh (2.12) 
T=0 


To consider heat losses, let us first have a look at the conductive losses for a pipe filled 
with a liquid at one temperature T'"#id. Here, the heat losses q'°® can be described 
with 

ge = kyndL Gree = pae) (2.13) 


using diameter d and length L of the pipe, ambient temperature T™”? and heat loss 
factor ky. This loss factor can be found in the pipe data sheet or be calculated based 
on material and construction. 


In a dynamic setup the temperature inside the pipe is not constant. Thus, to calculate 
the pipe outflow temperature T?’ as a function of the ambient temperature T™”? and 
past pipe inflow temperatures, the heat (loss) balance for one infinitesimal volume el- 
ement with length dx traveling through the pipe at position x needs to be calculated. 
Now let us combine (2.9) and (2.13) and adapt both for an infinitesimal volume ele- 
ment dx. With T(x) as the temperature of the infinitesimal volume element at position 
x and dT(x) as the change of this temperature, we get 


rdky dx (T(x) = re) = ritcpAT(x). (2.14) 
Solving this in-homogeneous differential equation for T(x) leads to 
ZEV, 
yo — Tamb ef (i = 2) e Pepa Tt (2.15) 


for the pipe outflow temperature T?"! depending on ambient temperature T""?, the 
travel time needed to pass the pipe tT, the past inflow temperature T;",, and pipe 
parameters like heat loss coefficient ky and diameter d. 


In a discrete time setting, the travel time T; can be calculated based on the number 
of time slots needed to pass the pipe. Using a simulation environment supporting 
continuous time equations like Modelica, (2.15) can be implemented directly in com- 
bination with 


Ut 
TH1— T= 1— 


2.16 
un (2.16) 


calculating the travel time 7; based on the velocity v; [SFT*05, Girl6]. 


Extensions of the node method 


An extension of the node method to a continuous time model is presented in [SFT*05, 
Girl6] which enables exact formulations in the model and simulation environment 
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Modelica. Oppelt presents a refinement of the discrete time node method in [Opp15], 
which tracks the propagation of temperature not only across one but across several 
pipes. Thus, averaging of temperatures occurs less often and results get more accu- 
rate. However, this refinement of the node method addresses simulation and is not 
straight forward to use in optimization of thermal grids. 


2.1.3 Storing thermal energy in a heating grid 


If the heating grid dynamics introduced in the previous section are used in a smart 
way, thermal energy can be stored in a heating grid. In Example 2.2, this storage effect 
is explained for a simple heating grid with one producer and one load. This example 
follows our publication [MHH19]. 


Example 2.2: 


We consider the most simple setup of a heating grid with one producer and one load 
connected with a supply and a return pipe as shown in Figure 2.4. A supply temperature 
profile at the producer is given with an increase of temperature in time step 3 and a decrease 
of temperature in time step 8. Perfect control of return temperature is assumed at the 
consumer leading to a constant return temperature. In average the supply temperature 
is slightly lower at the consumer due to thermal losses in the pipe. A constant thermal 
demand is assumed for the full time horizon. 


As the flow velocity in the pipe is limited, the temperature increase reaches the consumer 
with a transport delay. Hence, despite the temperature increase in time step 3 the mass 
flow first remains unchanged which leads to an immediate increase of thermal generation 
of the producer according to (2.9). As soon as the increased supply temperature reaches 
the consumer after two time steps in time step 5, the controller at the consumer reduces 
the mass flow to keep (2.9) in balance. As hydraulic dynamics are very fast (see Section 
2.1.2) the mass flow is reduced almost immediately at the producer, leading to a reduction 
of thermal energy output in time step 5. Thus, in time step 3 and time step 4, when the 
increase of temperature did not reach the consumer yet, the thermal energy produced at the 
producer exceeds the thermal demand of the consumer and the thermal losses. Hence, here 
the thermal dynamics of the heating grid show a behavior similar to charging an energy 
storage. 


Reducing the supply temperature at the producer in time step 8 leads to an inverse behavior. 
As the decrease of supply temperature is not seen at the consumer before time step 10, the 
mass flow stays constant until time step 10 and, thus, the thermal energy output of the 
producer is reduced proportionally to the supply temperature decrease. As soon as the 
supply temperature decrease reaches the consumer, its controller will increase the mass flow 
in order to supply the constant heat load. With this additional mass flow, the thermal energy 
output of the producer is increased back to its original level. In the period of time step 8 
and 9, the thermal energy output of the producer is below the thermal load plus thermal 
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losses. Thus, it shows a similar behavior as a discharge of an energy storage. Hence, the 
thermal dynamics of a heating grid show an energy-storage-like behavior where increases 


of supply temperature lead to charging of the storage and decreases of supply temperature 
lead to discharging of the storage. 
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Figure 2.4: Storage behavior of a heating grid with one CHP and one consumer. a) Supply temperatures 
at producer and load. b) Heat production of producer and heat consumption of consumer. 


Remark: As the mass flow is reduced in times of high supply temperature at the consumer, 
in reality the transport delay increases. Thus, the decrease of temperature in time step 8 


would arrive slightly later at the consumer. For high temperature differences this effect is 
more important than for small temperature changes. 
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2.1.4 Optimization of heating grids using the node method 


To utilize the storage effect explained above, several approaches in literature have 
been proposed. Unfortunately, the mass flow dependent time delay induced by the 
transport time in (2.15) as well as the bilinear terms in (2.9) lead to a non-convex 
problem. Thus, finding optimal solutions is very complex. 


A common approach used in many publications is to fix the mass flows or time delays 
prior to optimization, as this leads to a linear and thus convex problem. Some publica- 
tions assume directly that the mass flows in the heating grid are known upfront. This 
is the case when they are operated with a different control strategy without valves at 
the consumer side and, therefore, no changes in mass flow are induced by the con- 
sumers [LWW116, GWL*17, LWLL17, ZFZC18, ZZZW18]. Unfortunately, assuming 
that consumers do not control the flow means that no thermostats are used, which 
does not allow accurate local control of the indoor temperature. Hence, assuming 
constant mass flow can be a valid assumption in some areas, but for most heating 
grids worldwide this assumption is not valid. 


To overcome this issue, several researchers propose to use an iterative or sequen- 
tial solution approach. The authors of [BBR95] use a model with fixed time delays 
in iteration with a model with fixed supply temperatures to optimize operations of 
heating grids aiming for a fully dynamic optimization. Another iterative scheme is 
used in [SFT*05]. Here parameters of an optimization model with fixed time delays 
are iteratively updated using a simulation model. Robustness of the solution is in- 
creased by requiring an additional heat injection. A similar iterative approach is used 
in [Girl6é, GMBV17] where results of a detailed simulation update the parameters 
of a linear grid model. To get a linear model, fixed time delays from generation to 
consumers are assumed for optimization. In [DCM+19] the authors as well use an 
iterative optimization scheme where an optimization model with fixed delay times is 
updated using a simulation of the heating grid. In addition, the thermal inertia of 
buildings is considered as thermal storage opportunity . 


An intra-day optimization approach using non-linear optimization models is pro- 
posed in [LWS*16]. As intra-day approach it does not decide on generation unit com- 
mitment and does not consider forecast uncertainty. To cope with transport delays it 
uses an approach different from the papers mentioned before: Several complicating 
variables are identified for a pipe model based on the node method (including inte- 
ger time delays). Those complicating variables are fixed in a non-linear optimization 
model which is updated using a simulation model in every iteration. Similarly, the 
approach proposed in [Trö99] does not use fixed time delays. It relies on the prin- 
ciple ideas of the node method for deriving suitable models, but uses an iterative 
approach combining a static hydraulic, a static thermal and a dynamic thermal model 
to optimize the dynamics of a heating grid. 


14 2 State of the art 


2.1.5 Optimization of heating grids: other approaches 


Despite many publications using the node method to model the thermal dynamics 
of a heating grid, there are a some papers using different approaches. Starting from 
a general thermal energy balance a PID control approach considering the inherent 
storage capacity of a heating grid is presented in [BJPS11]. No optimization, but a 
parameter study is used to identify good times for charging or discharging the district 
heating grid’s storage capabilities. Lesko and Bujalski propose to model thermal 
dynamics of heating grids using a simple energy balance [LB17]. Here, the heating 
grid is assumed to be one single water volume with perfect mixing. Increasing or 
decreasing the supply temperature charges or discharges the thermal capacity of this 
volume. This modeling approach is compared with a model derived by the node 
method and parameterized for real-world data. It is shown that accuracy of both 
modeling approaches depends on a good parameterization [LB17]. 


In [Gro12], a rather different approach is used: The heat storage capabilities of a 
radial heating grid are approximated with a linear regression model which is trained 
based on a high number of heating grid simulations. Those grid simulations require 
a detailed physical model of the grid, whereas the resulting linear regression model is 
a lightweight model linking the thermal energy stored inside the grid with past and 
current supply temperatures and the current load factor of the grid. Hence, it enables 
fast optimization runs as the computational effort is moved from optimization to 
model parameterization. 


In [VD15], a detailed simulation model built in MATLAB/Simulink derives heat de- 
mands and mass flows in a meshed heating grid. Those results are used as inputs 
for an optimization of the generation sites. This work is extended in [VTD17] with 
a hybrid-evolutionary algorithm combining a MILP at lower level and a genetic al- 
gorithm at upper level to minimize total operating cost of a meshed heating grid. 
However, thermal dynamics of heating grids are only considered in steady state in 
the simulation model run by the genetic algorithm. 


A decomposition of the mixed-integer non-linear problem (MINLP)? of district heat- 
ing optimization into a unit commitment and an economic dispatch problem is pro- 
posed in [SLM*17]. The unit commitment problem is formulated as mixed inte- 
ger quadratic constraint program without consideration of the heating grid dynam- 
ics leading to a problem similar the one presented in Section 2.1.1 except having a 
quadratic instead of a linear cost function. The economic dispatch problem is a non- 
linear program (NLP)? which is formulated in the object-oriented modeling language 


Like MILPs, MINLPs are optimization problems consisting of real and integer variables. However, 
for MINLPs non-linear objective functions as well as non-linear equality and inequality constraints are 
possible. As MINLPs are hard to solve, off-the-shelf solvers only exist for special types of non-linearity. 
In difference to MINLPs, NLPs are optimization problems only consisting of real variables. Off-the-shelf 
solvers finding local optima are available for NLPs. 
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Modelica and considers full dynamics of the heating grid. A local solver is used to 
solve this NLP. 


As methods for optimization of district heating grids are usually computationally 
expensive, there have been efforts to speed up simulation and optimization of heating 
grids. A Danish and a German method for aggregation of the grid topology trying to 
loose as little accuracy as possible have been proposed. A comparison of the two is 
published in [LBW04, WBS05]. 


Considering different forms of energy beyond district heating, Geidl and Andersson 
propose a joint optimal power flow for electricity, heat and gas networks [GA07]. 
To enable applicability to large scale systems, a decomposition scheme for the afore- 
mentioned approach is presented in [MAAFFH14]. However, for both publications 
heating grid dynamics are not in focus and variable transport delays in temperature 
propagation are not considered. In [LM17, ZFZC18], a combination of electricity, 
gas and heat networks is considered as well, whereas the building thermal inertia is 
used as additional thermal energy storage besides the heating grid thermal inertia in 
[LWLL17, GWL+17]. 


A completely different storage concept is proposed in [DMOK19]. Here, instead of 
using the thermal dynamics induced by variations of supply temperature, the flow 
direction is reversed in order to store thermal energy inside a heating grid. This 
approach for thermal storage is only applicable if large thermal generation exists at 
both sides of a pipe used for storage. Hence, it implies strong requirements on system 
setup and grid topology and, therefore, is not considered further in this thesis. 


2.2 Global optimization of non-convex problems 


Global optimization of non-convex problems is a complex topic studied for many 
decades [HT96]. A common, well-established approach proving global optimality of 
a solution, is to transform the non-convex problem into a primal and a dual prob- 
lem which both are easier to solve than the original non-convex problem [FV93]. The 
primal problem leads to a feasible solution of the original problem, for example by 
searching for a solution only inside a reduced solution space. For a minimization 
problem, this solution of the primal problem represents an upper bound of the origi- 
nal problem, as better solutions might be possible in the original solution space. The 
dual problem is an outer-approximation of the original problem. Very often it is 
chosen to be convex in order to ensure global optimality of its solution. It increases 
the original solution space such that all feasible solutions of the original problem are 
part of the solution space of the dual problem, but additional solutions are possi- 
ble. Hence, for minimization problems, the solution of the dual problem represents 
a lower bound of the original problem. Solutions of the original problem cannot be 
better, but the dual solution might not be a feasible solution of the original problem. 
To solve the original non-convex problem, the primal and dual problem are solved 
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in multiple iterations. If primal and dual solution are equal, the globally optimal 
solution is reached. If they are not equal, the primal solution is a feasible solution 
of the original problem and the dual solution indicates the gap to global optimality. 
If in every iteration the exactness of the dual and primal solutions is improved, this 
gap is reduced. The iterations usually stop, when a specific numeric precision of the 
solution is reached (e.g. a gap smaller than 1 %) [FV93]. 


For some non-convex problem classes like MILPs, solvers using such a primal-dual 
algorithm reaching proven global optimality are available as open source [BENO5, 
FRV*+18] or commercial tools [Gur16]. 


In the following, only the problem classes relevant for modeling and optimization 
of district heating grids are discussed where no off-the-shelf solvers exist. These are 
bilinear terms like in (2.9) or (2.10) and a variable-dependent time delay as in (2.15). 
Several approaches exist for global optimization of bilinear problems as discussed in 
Section 2.2.1 and Section 2.2.2. Necessary conditions for a local optimum of variable- 
dependent time delays are defined in [CGCP16] and a solution algorithm finding such 
local optima is proposed in [CPB17]. However, for problems with variable-dependent 
time delays, no approach proving global optimality of its solution has been presented 
yet. Global optimization solvers for NLPs or MINLPs like BARON [Sah96] do also 
not support solving variable-dependent time delays. 


2.2.1 McCormick envelopes 


McCormick envelopes are commonly used as outer-approximation or dual problem 
for global optimization of bilinear terms with bounded multipliers, because they are 
the tightest convex hull of a bilinear term with bound variables [McC76]. Figure 2.5 
shows the product of two variables, each bounded by 1 and 5, with under-estimating 
McCormick envelopes. It can easily be seen that combining the planes in the two 
sub-figures provides a tight convex hull of the bilinear term. In combination with a 
branching strategy or a piece-wise linear approximation, McCormick envelopes can 
be used to build a global optimization scheme for bilinear terms. 


Therefore, McCormick envelopes are employed in several approaches for global opti- 
mization of bilinear terms. In [AKF83], a branch and bound solution algorithm using 
McCormick envelopes as over- and under-estimators of a bilinear function is pre- 
sented. A piece-wise McCormick relaxation is used in [4ACZ*17] for handling bilin- 
ear terms in a scheduling problem for operation of crude oil terminals. A piece-wise 
linear relaxation with McCormick envelopes using bivariate partitioning is presented 
in [HK10] to solve bilinear programs. 
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Figure 2.5: Surface plots of product w = a-b with a,b € [1,5] and 
planes of the McCormick under-estimators 


2.2.2 Multiparametric disaggregation 


In this thesis, multiparametric disaggregation [TCM12] is used for global optimiza- 
tion of bilinear terms as it promises faster solution times than a piece-wise linear 
approximation using McCormick envelopes [CT13]. 


The concept of multiparametric disaggregation is the discretization of one of the vari- 
ables in a product w;; = xi: xj. Like with floating point numbers in computing, every 
digit of this variable is split into multiple binary variables which encode the value of 
this digit. Here, variable x; is split into binary variables z;,; with I representing the 
position of the digit and k taking values from 0 to 9 representing the value of this 
digit. Hence, as every digit can take only one value, z;x; is non-zero only for one k for 
all! and j [TCM12]. 


P 9 
a 2.2 10'-k- Ziki + AX; (2.17) 
i=p k=0 
9 
Yau-ı AMezeziep (2.18) 
k=0 
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As | defines the position of the digit, p and P define the precision of approximation. 
If p approaches negative infinity and P is large enough, (2.17) becomes exact. Ax; is a 
slack variable bounded by 0 < Ax; < 10? which transforms the product Wij according 
to [TCM12] into 


P 9 
KEE, L zw -k ĉiii + Awjj. (2.19) 
=p ae 


Here ijn = Xi Zjkl which is enforced with 


9 
k=0 


xP + zu S Rijn <x zr {Yke ZOOSK <9}, {VIE Z|p<1<P} (221) 


1 


where xt and a are the lower and upper bounds of variable x; [TCM12]. 


The slack variable Aw;; is linked to the slack variable Ax; by Aw; = xi: Ax; and 
accordingly [TCM12] 
xh Ax; < A wy < x} - Axy. (2.22) 


Due to the slack variable Aw;; the solution space of (2.19) includes the original solution 
of the product w;; = x; x; and, thus, it is an outer approximation [TCM12]. Hence, 
this yields a MILP which is a dual (or lower bound) problem of the original problem. 
The iterative scheme shown in Figure 2.6 is used to solve the bilinear problem to 
global optimality, using the non-linear model as primal (or upper bound) problem 
being solved with a local, non-linear solver [TCM12]. 


For ease of understanding, multiparametric disaggregation is introduced using a base 
of 10 (k = 0 to 9) in this section. Choosing a different base is possible potentially 
leading to faster solution times [TCM12]. 


2.3 Pipeline scheduling: modeling of 
flow dynamics in other areas 


In addition to the optimization approaches presented in Sections 2.1.4 and 2.1.5, opti- 
mization models for related problems can inspire innovative solutions for heating grid 
optimization. In this section, pipeline scheduling is introduced as it is a well stud- 
ied field of research with many similarities to heating grid optimization. In pipeline 
scheduling, liquid refined petroleum products need to be transported via a pipeline 
system from refineries to depots. There are multiple products, but only a limited 
number of pipes connecting refineries and depots. 


In Figure 2.7 a pipeline scheduling example presented in [MC17] is shown. There are 
three pipe segments in this example. Pipe segment S1 connects the input node with 
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Figure 2.6: Iterative solution approach of multiparametric disaggregation 


refinery R1 to the output node with depot D1. Segment 52 connects the output node 
with Depot D1 to the dual purpose node DP1 with input by refinery R2 and output 
to depot D2. Pipe segment S3 connects dual purpose node DP1 to the output node 
with depot D3. There are five different products P1 to P5 which can be produced at 
refineries R1 and R2. Depots D1, D2 and D3 each require a certain amount of these 
products. The goal of optimal pipeline scheduling is to find an optimal product se- 
quence transported in the pipes allowing to achieve all demands at minimum cost 
or in minimum time. Pipeline scheduling for refinery products is having many simi- 
larities with optimization of heating grids: There are nodes feeding into the pipeline 
system and nodes consuming from the pipeline system. The dynamics of propaga- 
tion of a liquid in a pipe are similar, too. The liquid propagates with the mass flow 
and there is a transport time from entering to leaving the pipe. However, in pipeline 
scheduling for refinery products there is a discrete set of different products which is 
transported via the pipes, whereas in heating grids the supply temperature is a real 
variable which can be chosen freely within its limits. 


In [RP04], a discrete volume representation of the pipeline in combination with a 
discrete time grid for the scheduling horizon is introduced to optimize the pipeline 
scheduling problem for refinery products. However, using a continuous time for- 
mulation of flow dynamics, [CC04] allows a MILP formulation. Today’s models for 
pipeline scheduling are based on this concept and can handle reversible mass flows 
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Figure 2.7: Pipeline distribution system for refined petroleum products after [MC17]. 


and different network topologies using a continuous time formulation as well as a 
continuous representation of volumes [Cas17]. A major difference of this product- 
centric model formulation to other model concepts like [RP04] is that the variable flow 
rates are considered only implicitly without being an explicit optimization variable. 
This allows a very efficient MILP continuous time formulation [CGZ18] avoiding non- 
linearities in the optimization. However, it prevents an accurate modeling of pumping 
cost [CCMC15]. 


2.4 Objective of this thesis 


CHPs need to operate more and more flexible to react to volatile renewable genera- 
tion, volatile demand and volatile energy prices which requires thermal storage. The 
thermal dynamics of a heating grid allow to store thermal energy. Thus, they enable 
a flexible operation of CHPs without requiring an investment in dedicated thermal 
storage tanks. To leverage this potential, the thermal dynamics of a heating grid need 
to be considered in operations planning. Model formulations for these dynamics 
usually include bilinear terms and a variable-dependent time delay representing the 
temperature propagation in a pipe. Thus, they result in non-convex problems and the 
approaches proposed to solve these problems do not prove the global optimality of 
their solution as shown in Section 2.1.4 and Section 2.1.5. In the literature on global 
optimization, there are many approaches for bilinear problems (see Section 2.2.1 and 
Section 2.2.2). However, there are no approaches for global optimization of problems 
with variable-dependent time delays. 


Thus, the approaches for heating grid optimization presented in the literature do not 
offer a guarantee to reach the global optimum and the literature on global optimiza- 
tion does not offer appropriate solution algorithms either. Nevertheless, it is essential 
to find a solution as close as possible to the global optimum to operate a heating 
grid as sustainable as possible and to leverage the full potential of using the thermal 
dynamics of a heating grid as energy storage. Accordingly, this thesis aims at find- 
ing the globally optimal solution for scheduling of CHPs considering heating grid 
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dynamics. 


In addition, current optimization approaches considering the thermal dynamics of 
heating grids are only very rarely applied in real world installations. Thus, this the- 
sis aims as well at finding real world applicable optimization models which enable 
fast and reliable solutions. A model training using historic measurements would be 
ideal to enable easy brown field installations without a major model parameterization 
effort. 


In summary, the following goals in modeling and optimization of heating grids will 
be addressed in this thesis: 


e Find a globally optimal solution for optimization of heating grids allowing to lever- 
age the full potential of considering the thermal dynamics in operations planning 


e Find fast and reliable approaches for optimization of heating grids with accurate 
solutions and easy model parameterization enabling real world applications 


3 Global optimization with 
multiparametric delay modeling 


Existing approaches as discussed in Section 2.1 do not prove global optimality of their 
solutions as they simplify the non-convex heating grid dynamics in their optimization 
models. Hence, it is not possible to judge the quality of their solution as it is unknown 
if a better solution exists. Non-convexities in heating grid dynamics are bilinear terms 
in (2.8) and (2.9) as well as a variable-dependent time delay in (2.15). Global optimiza- 
tion of bilinear terms is a well studied field as discussed in Section 2.2. However, there 
are no methods for global optimization of problems with variable-dependent time de- 
lays. In the following, “multiparametric delay modeling” is presented, that is able 
to deal with this problem using a primal-dual algorithm as discussed in Section 2.2. 
Starting with the dual problem in Section 3.1, this chapter outlines two methods to 
formulate the primal problem in Section 3.2 and proposes the resulting primal-dual 
algorithm in Section 3.3. 


“Multiparametric delay modeling” is the first approach for global optimization of sys- 
tems with variable-dependent time delays. It enables to globally optimize the fully 
dynamic problem of optimal operation of heating grids if combined with a global op- 
timization approach for bilinear terms like multiparametric disaggregation presented 
in Section 2.2.2. The basic idea of “multiparametric delay modeling” was first pub- 
lished in [MH19]. More details are published in [MLH19] with a discussion of the 
primal and dual problems and a presentation of proofs of convergence to the global 
optimum. This chapter presents the enhanced model formulation as in [MLH19] as it 
improved significantly solution quality and solution times. 


The optimization model for global optimization of heating grids is based on the unit 
commitment problem without thermal dynamics (Section 2.1.1) where the heat bal- 
ance equation (2.5) is replaced with the dynamics of the node method (Section 2.1.2). 
The following sections discuss how to model the non-convexities of this problem for- 
mulation in the primal and the dual problem. 


3.1 Dual problem 


As a first step to a primal-dual global optimization algorithm, this section introduces 
the outer approximation or dual problem of variable-dependent delay times based on 
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the node method (Section 2.1.2) as in [MLH19]. The basic idea of this outer approx- 
imation of the pipe outflow temperature with “multiparametric delay modeling” is 
explained in Example 3.1. 


Example 3.1: 


As in Example 2.1, Figure 3.1 shows a temperature distribution inside a pipe at beginning 
of time slot t. The pipe of length L is filled with volumes at temperature levels of past 
inflow temperatures pee T and Tha Those volumes span a width of Ax;_1, Ax+—2 
and Ax;_3 being the propagation of water inside the pipe in the previous time slots t — 1, 
t—2andt—3. 


Figure 3.1: Pipe filled with different temperatures at beginning of time slot t. The red-dashed volume 
ji leaves the pipe in time slot t. The blue arrow indicates the flow direction (left to right). 


The red-dashed volume jı in Figure 3.1 will leave the pipe in time slot t. Thus, its average 
temperature is the average pipe outflow temperature in this time slot. The temperature of 
this red-dashed volume jı is a mixture of volumes being at temperatures T}", and Ti" 4. 
Hence, the average outflow temperature T?"' in time slot t needs to be between those two 
temperatures. Thus, the following inequality leads to an outer approximation of the outflow 
temperature Te“ in time slot t. 


min (Ty, Tis) < Te < max (Te Tits ) (3.1) 


To refine this outer approximation, the volume leaving the pipe in time slot t is divided into 
two sub-volumes jı and ja of identical size in Figure 3.2. 


in in in 
Ti- Ta Tı-3 


Figure 3.2: Pipe filled with different temperatures at beginning of time slot t. The volumes jı (red) and 
j2 (green) leave the pipe in time slot t. The blue arrow indicates the flow direction (left to 
right). 
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Now the outflow temperature of sub-volume j; (red) Tee is known exactly as it only con- 


sists of water at temperature Ti". The outflow temperature of sub-volume jz (green) Tee 
can be estimated using a similar outer approximation as before for volume jı in Figure 3.1: 


Tpit = Tia, (3.2) 


: in in out in in 
min ( t-2, dt 3) < Tih < max ( t-27 i). (3.3) 


To get the pipe outflow temperature Te“' in time slot t assuming perfect mixing of sub- 
volumes jı and ja, the average of the temperatures of the sub-volumes needs to be calculated 
with 


1 1 


ees ae (3.4) 
All in all, for the case shown in Figure 3.2 we get 
1 in > in in out 1 in 1 in in 
2 t 5 min ( t-2 Kt 3) <T" < alt 2t > max ( t-2 ey (3.5) 


as outer approximation for the pipe outflow temperature T?"! in time slot t. 

In (3.5) the temperature of only half of the volume leaving the pipe in time slot t is estimated. 
Hence, (3.5) is a tighter outer approximation than the outer approximation in (3.1) and 
adding a sub-volume increased tightness of the outer approximation. 


Example 3.1 shows a way to formulate outer approximations for the outflow tem- 
perature of a pipe including a possibility to refine this outer approximation by intro- 
ducing sub-volumes. Thus, if there is a way to formulate this approach in general 
as a solvable problem, it can be used as dual in a primal-dual global optimization 
scheme. The following presents a general MILP formulation of this concept as in 
[MH19, MLH19]. 


The non-negative real variable Ax; represents the distance traveled by the transport 
medium inside the pipe in one time slot t. It is calculated using (2.11): 


er eh (3.6) 


or (4)" 


Based on this distance traveled Ax;, we introduce the binary variable b, ; ; representing 
a discrete time delay. This binary variable bq ; + shall equal 1 if and only if sub-volume 
j leaving the pipe in time slot t entered the pipe T time slots ago. Hence, for every 
time slot t and every sub-volume j exactly one time delay T exists and all other bz j, 
have to be zero, which can be enforced by 


Ax; = At -v = 


t 
L br it = Vt € S; Vj € Sj. (3.7) 
T=0 
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Here, Sj denotes the set of sub-volumes j. This set is limited by the number of sub- 
volumes np: Sj € [1, 1p]. 


To ensure that the correct time delay T is selected, the following inequalities are intro- 
duced. They compare the sum of past distances traveled within the pipe (Axı) with 
the length of the pipe L adjusted by the sizes of the j sub-volumes. Using a sufficiently 
large parameter M* allows to set the correct discrete delay variable bq j; to 1 for every 
sub-volume j in every time slot t with 


i—1 t—1 
R= L Au< E Axut (1-b)  M*  YEESpYjES YTES, (88) 
P t=t—-T 


—1 t—1 
L- An L Are- (1-bea) M* VEESLViES, VTE Sr. (9) 
np t=t—T+1 


The set S; contains all possible (discrete) time delays. 


Next, we need to find suitable outer approximations of the outflow temperatures of 
these sub-volumes TH. Here several different cases need to be considered depending 
on the number of possible temperatures within a sub-volume. 

First, if the temperature inside the sub-volume is constant, as for sub-volume jı in 
Figure 3.2, the outflow temperature of the sub-volume Ten is known exactly. This 
is the case when the discrete time delays b,j; of two neighboring sub-volumes are 
equal. Thus, with TU being the upper bound of the temperature we can assign the 
correct past inflow temperature Tj”, to the outflow temperature of the sub-volume 


Tp using 


TH < T+ (2— ihre) TU Vt E Sp Yj € Sj}, Vt € Sr (3.10) 
Tr STi (2b bypass) TY VE E Se Vj € Sj, YT € Sr. (8.11) 


Second, two different temperature levels are possible, if the discrete time delays bq j, 
of two neighboring sub-volumes j and j + 1 are different by 1 (br; and b-_1,j+1,: are 
non-zero). Thus, an outer approximation of the outflow temperature as proposed in 
Example 3.1 can be used leading to 


Tpit < max (ms nu) +(2- br jt — br-1,j41,4) TY 
Vt © Si, Vj = Sj, YT pze Sr (3.12) 


ty men (Tin Ti u 
Tj 2 min (ies 2 = (2 — by jt — brij) T 


Vt = St Vj = Sj,VT Se Sr. (3.13) 


3.1 Dual problem 27 


A linear formulation allows an efficient solution of a MILP problem. Thus, the func- 
tions min(a,b) and max(a,b) need to be represented in linear equations. For this, we 
note that for any real, non-negative variables a € R>o and b € R>ọ the difference 
a — b can be decomposed into two real, non-negative auxiliary variables ö* € Rso 
and ö” € R>o as follows: 

a—b= ôt — ô (3.14) 


Using these auxiliary variables we can enforce c < max (a,b) with 

e<b+öt. (3.15) 
Similarly for c > min (a,b) we get 

c>b-6”. (3.16) 
Third, if discrete delay times br jt of two neighboring sub-volumes j and j + 1 have a 
distance of 2 (b;,;; and b;_2 j41, are selected), the sub-volume outflow temperature is 
a mix of three temperatures. Hence, the sub-volume outflow temperature needs to be 


between the minimum and the maximum of those three temperatures. In line with 
the outer approximation above we get 


Tp < max (m z Ti ay I. : 2) + (2 = bzr jt — br-2j+1,) . TU 
Vt © S, Vj kaes Sj, YT ja Sr (3.17) 


a ee | a 
Tptt > min (Th, Th gt Tet ega) — (2— br je — be ajaae) T 


Vt St, Vj = Sj, YT S Sr. (3.18) 


The formulation from (3.14) to (3.16) for min and max functions can be extended to 
calculate the minimum or maximum of three variables using 


min (a,b,c) = min (min (a,b), c), 


max (a,b,c) = max (max (a,b),c). 


The implementation of cases with more than three temperature levels in a sub-volume 
is not required. Indeed, by increasing the number of sub-volumes n, we always obtain 
a situation with at most two temperatures per sub-volume as proven in Theorem 3.1. 
In early iterations of the optimization algorithm sub-volumes are still large. Thus, 
they easily span more than two or three temperature levels. To increase tightness 
of these early approximations, it is recommended to limit the sub-volume outflow 


temperature (tT < Te < su) for all sub-volumes j and time slots t. In fact, the 
outer approximation for three temperature levels in (3.17) and (3.18) is theoretically 
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not required to reach the global optimum as with a high number of sub-volumes this 
situation does not occur. However, it increases tightness of the formulation in early 
iterations, too, and this tighter formulation reduces solution time. 


Theorem 3.1 (Decreasing number of temperatures per sub-volume) 


There is an n, € IN such that for all np > nj, there are at most two temperatures inside 
any sub-volume j € Sj. 


Proof: 

To prove Theorem 3.1, first note that the inflow temperature T?” into a pipe only 
changes at the beginning of a time slot t having a discrete time model. Thus, the 
distance between two neighboring temperatures T/", and TÌ” in a pipe is the progress 
of the transport medium within one time slot Ax; which can be calculated according 
to (3.6). 


With the proposed approach, the volume leaving the pipe in time slot t is divided into 
np equal-sized sub-volumes j. As the volume leaving the pipe in time slot t spans Ax;, 
the np, sub-volumes j leaving the pipe in time slot t span Ax; ;: 


Ax f 
Agee m Vj E SYE E Si. (3.19) 


If the number of sub-volumes n, is increased, the span of one sub-volume AxXjt de- 
creases. Without any limitation on Ax; we arrive at 

f _ Axı : 

lim Axj;= lim —— =0 Vj E Sj, VE E Se. (3.20) 


No Np—oo Np 


Thus, there is always a number of sub-volumes n} big enough to guarantee 


Axjp< Axi VEE Se {Vit € Sitt < t}. (3.21) 


If more than two temperature levels are contained within a sub-volume j at time t, 
there exist two temperatures in the past (Vtt € S,|tt < t) having a smaller distance 
Axy than the size of the sub-volume Ax; . However, if the number of sub-volumes np 
is increased above n}, inequality (3.21) holds for all sub-volumes j and time slots t and, 
hence, more than two temperatures within any sub-volume j are not possible. 


To calculate the pipe outflow temperature T?” in time slot t, the average of the tem- 
peratures Te of the np equal-sized sub-volumes j leaving the pipe in time slot t needs 
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to be calculated like in Example 3.1. As outer approximation of the pipe outflow tem- 
perature T?"! we get 


de 
out = out 
Tot > I = Te veS 3.23 
pe np jZ bj San (3.23) 


As motivated in Example 3.1 and proven in Theorem 3.2, an increase of the num- 
ber of sub-volumes n, leads to an improved precision of this outer approximation. 
Thus, the outer approximation for variable-dependent time delays not considering 
thermal losses presented in (3.10) to (3.23) can serve as a dual in a primal-dual global 
optimization scheme as introduced in Section 2.2. 


Theorem 3.2 (Increasing tightness of outer approximation) 


If the number of sub-volumes approaches infinity (np — 00), the approximation of the 
pipe outflow temperature presented in this section approaches its real value Te". 


Proof: 

To prove Theorem 3.2, we note that the pipe inflow temperature T?” only changes at 
the beginning of time slots t as a discrete time modeling approach is used. Hence, 
at every time slot t, the number of past inflow temperatures is limited to the num- 
ber of past time slots t — 1. Introducing a variable for the number of past inflow 
temperatures nr, in time slot t we can state accordingly 


nT <t<o Vt E Si. (3.24) 


In reality, the number of past inflow temperatures nr, is even more restricted, as 
there are physical limitations on mass flows and velocity of the transport medium. 
The number of possible temperatures inside a sub-volume leaving the pipe is bound 
by the same number nr +, as only past inflow temperatures can be inside a sub-volume 
leaving the pipe. 


Increasing the number of sub-volumes n, we will always reach a situation where 
Np > NTs (3.25) 


Hence, there are np — nr,; more sub-volumes than possible temperatures and for np — 
nr, sub-volumes there is at most one temperature possible. Thus, for these np — nr, 
sub-volumes the sub-volume outflow temperature Tp is known exactly. Having a 
limited nr, as shown in (3.24) for large np we get 


Np NT > NT 4. (3.26) 
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Recalling (3.22) and (3.23), the overall outflow temperature T?"! in time slot t is the 
average of the temperatures of sub-volumes Tp. 


i 
Tet = z = Tor VEES: (3.27) 
j=l 


Increasing the number of sub-volumes np, the at least np — nr,; sub-volumes with 
exactly known outflow temperatures ur get increasingly more weight in (3.27) than 


the at most nr; sub-volumes with approximated outflow temperatures Ten. Thus, 
the outer approximation is getting more tight and exact with increasing np and 


n 


1 
lim — )° Tp = Tet tE Si 2 
ns F L T h Vt E S; (3.28) 


So far, (3.10) to (3.13) and (3.17) to (3.18) approximate the outflow temperature of one 
sub volume without consideration of thermal losses. To integrate thermal losses, these 
equations need to be combined with (2.15). Thereby, (3.10) becomes 


; -V At 
tS Te eee = ye") eT + (2— beje — De jst) TY 
—_ E 


.— Fin 
=T 


Vj SS Sj, Vt < Si, YT La Sr. (3.29) 


As in (2.15), ky, cp, p and d are parameters of pipe and transport medium, At is 
the length of one time slot and T is the discretized time delay coming with bz j. 
For the remaining (3.11) to (3.13) and (3.17) to (3.18) approximating the sub-volume 
outflow temperatures, heat losses can be integrated accordingly with Ti”, as defined 
in (3.29): 


TE > Tit, (2-— brjt— de jyt) TY, VEE SE ViES; YTES: (3.30) 


t Ti Fi u 
Tij < max (Tiy ag) + (2 ~~ br jt = br1,j+1,) -T 


Vt © St Vj h Sj, YT < Sr (3.31) 


ty in (Fin Fi u 
Tj 2 min (Tes, as) = (2 = by jt -br-1,j+1,)T 


Vt =, Si Vj = Sj,VT = Sr (3.32) 
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Tp < max (Ti u Hr To. 2) + (2 _ br jt _ br_2;41,) . Tu 
Vt © S, Vj S Sj, VT | san Sr (3.33) 


out : Tin Tin Pin u 
Te > min (Dr Theat Te rag) — (2— br jt — ben) T 


Vt © St, Vj paa Sj, YT S Sr. (3.34) 


The min and max functions are implemented using (3.14) to (3.16). Consideration of 
heat losses does not change (3.22) and (3.23), the combination of sub-volume outflow 
temperatures to the pipe outflow temperature T?“. 


Combining (3.6) to (3.9), (3.22) and (3.23) as well as (3.29) to (3.34) we get an outer- 
approximation of the pipe outflow temperature as MILP problem using the basic 
concepts of the node method (Section 2.1.2). Thus, they can replace the equation for 
temperature propagation and loss (2.15) in the dual problem of a global optimization 
scheme for heating grids. To complete the MILP formulation of the dual problem, the 
bilinear terms in the node method required for temperature mixing in nodes (2.8) and 
change of inner energy at producers and loads (2.9) are modeled using multipara- 
metric disaggregation introduced in Section 2.2.2 and the unit commitment equations 
introduced in Section 2.1.1 are added. The resulting cost minimization problem of 
heating grid operations is summarized in Model 3.1. 


Remark: If the overall cost is minimized or profit or energy efficiency are maximized, 
the MILP solver will favor to produce less thermal losses. Thus, for heating grids the 
optimizer will maximize the pipe outflow temperature to avoid thermal losses and, 
hence, (3.22), (3.29), (3.31) and (3.33) are sufficient. If the approach is adapted for 
cooling grids only (3.23), (3.30), (3.32) and (3.34) need to be used as the optimizer will 
try to reduce the pipe outflow temperature to avoid thermal losses. (3.6) to (3.9), (3.22) 
and (3.23) are needed for both cases, heating and cooling grids. 


Model 3.1: Dual problem of multiparametric delay modeling 
Objective function 


3 . _el,b b ; 
min ), (prices! W. pi” — pricef!#ell. pit") 
tes; 


t el,var heat,var 
KE (ce Men > Pin +C; 44) 
tESt iESg 


Generation and load constraints 


nig PPS pis Sie PY Vi € Se, WE € 5; 
nig QF < git < nip QF Vi € Sg, Vt € St 
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b 
pre = D Dit + pr” Fr pel VEE S; 
i€ Se 
gis = opt, (Te — Tre) Vi € Sp U Sa VEE Si 


The constraints for the feasible regions of the CHPs need to be added matching 
their type (e.g. (2.3) for CHPs with extraction condensing turbine). 


Node balances 


L Mit = L Mit VEE St 


icin flow icout flow 


Lieinflow E rinis) 


Toit = - VES 
í Licinflow Mit i 
Assignment of binary delay 
At 
d 
pr (4) 
t 
insel VreSs,Vjes; 
T=0 
j Sai t—1 
L- j Axt < L AX + (1 ar br jt) - M* Vt € St Vj = Sj, YT Sr 
p H=t-T 


—1 t—1 
L-I—_Au> % Am-(l-bu):M Vees,Vjes,vres, 
np tt=t—T+1 


Pipe outflow temperature calculation 


Aky 
; — —IAtT 
He TUF (1 = roe) e Pept + (2— brjt— br jyt): TY 
O= 
=i, 


Vj Sj, Vt Si, YT Sr 


t Ti Fi u 
TO < max (Tr, Tr) +2 bejt — bern) T 


Vt € Si, Vj Sj, YT Sr 
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Ter < max (7 = Te SP TE | 2) + (2 — br jt = br-2,j+1,) . TU 
Vt © Si, Vj = Sj, YT Lsm Sr 


1 
Tu pout WEES 
en my H i 


All multiplications and divisions of variables are modeled using multiparametric 
disaggregation (Section 2.2.2). 


Function c < max (a,b) is enforced using the following equations. 


a-b=67-6 
c<b+6r 


3.2 Primal problems 


In addition to the dual problem presented in Section 3.1 a model for the primal prob- 
lem is required to complete a primal-dual global optimization scheme. Because of the 
variable-dependent time delay, it is unfortunately not possible to directly implement 
the node method as presented in Section 2.1.2 for non-linear optimization solvers 
like IPOPT [WB06]. Thus, to get a fast and easy to solve NLP as primal problem as 
in [TCM12], in the following two different approaches are introduced to model the 
variable-dependent delay. Both are first published in [MLH19] and adapt prior work 
[vFR*17, BHW17]. The approach presented in Section 3.2.1 uses first order differ- 
ential equations to model the transport delay. Section 3.2.2 presents a refinement of 
this approach combining a static time delay calculated using the node method with a 
variable delay model based on first-order differential equations. This leads to a sec- 
ond formulation for the primal problem. The bilinear terms remain for the primal 
problems presented in the following as they can be solved with local NLP solvers. 


3.2.1 Finite-volumes pipe formulation 


One common way to approximate dynamic time delays is to combine a series of first- 
order differential equations (PT1 elements) [vVFR‘17]. This approximation is getting 
more exact if the number of first-order differential equations increases. Applying 
this approach to the mass-flow dependent transport delay of a pipe, each first-order 
differential equation represents the energy balance of a volume element in the pipe. 
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GUS aera? 


0 n L 


Figure 3.3: Scheme of a pipe of length L split into n finite volumes i. 


As shown in Figure 3.3 the pipe is split into a finite number of n volume elements i of 
equal size. Assuming perfect mixing, the water inside each of these volume elements i 
is at only one average temperature T;’?. Hence, recalling (2.13) the thermal loss Cia 
of one volume element i in a time slot t becomes 


loss _ 
it 


kyrak (7 2) Vi € S YEE St. (3.35) 


As for (2.13), ky, d and L are parameters of the pipe and T“"? is the ambient tem- 
perature. The set S; denotes all volume elements. The transport of thermal energy 
Qe from one volume element i — 1 to the next volume element i is based on 
convection at mass flow m; with the transport medium having a specific heat capacity 
Cp and, thus, can be described with 


On — corte (Tis — TE) Wi € Si, YEE St. (3.36) 
The thermal energy balance for one volume element i is 
Qit =, Ore: = gig Vi € S, Yt € Si. (3.37) 


Combining the above equations with (2.9) for the inner energy change Q;; we get the 
energy balance of one volume element i with mass m = n/(o7(d/2)*L): 


kyrdL 
cpm: (Tifa — TE) = cprin (TH, — TH) — FE (Tg - r) 
Vi € Si, Vt E€ Sten, (3.38) 


With some minor reformulations, the resulting first-order difference equation is 


n kyrdL 
WaT + g (epr (TP, = Te) - AE (mp — rem) ) 
cpp (5) L 


N 
a 
NS 


Vi € Si, VEE Stem. (3.39) 


For a flow direction from first volume element i = 1 to last volume element i = n 
the outflow temperature of the pipe T?"' is the temperature T,, of the last volume 
element n. The inflow temperature of the pipe T!” replaces the temperature of the 
previous volume element T}, in the energy balance equation (3.39) for the first 
volume element i = 1. 
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The resulting model formulation is shown in Model 3.2. It allows to find fast and 
reliable solutions using local solvers for NLPs like IPOPT [WB06]. But as discussed in 
[Tr699] introducing a finite number of volume elements, each at only one temperature 
level, leads to increased averaging of the temperatures and, thus, smoothing of the 
delayed variable. Thus, it is less exact than the node method. This effect decreases 
with increasing number of volume elements n [Trö99]. However, a higher number of 
volume elements n increases problem complexity and solution time. 


Model 3.2: Primal problem using finite-volumes pipe model 
Objective function 


à . el,b b ; 
min ), (prices! YY — price!" | pi") 
tes; 


+E E (Gm + CPP pu cg.) 


tes; i€S¢ 


Generation and load constraints 


nit: PE < pit < nis PH Vi € Sp, Yt E Sı 
nit: QF < qin < nig: QU Vi € Sp VEE S; 
b 
paemand _ 3 Dit + Pi uy Spt VES: 
IE Sg 
Git = Cpi, (ria ER Vi € SgUSq,Vt € St 


The constraints for the feasible regions of the CHPs need to be added matching 
their type (e.g. (2.3) for CHPs with extraction condensing turbine). 


Node balances 
L Mit = 3 Mit Vt € St 
icin flow icout flow 


Licinftom (TH it nie) 


Lieinflow Mi, 
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Temperature propagation in pipes 
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ky7tdL 
Tir = Tir + z= (com (Tiit Ti) Z (1? 7) 
L 


Vi € Sizi, VE € Sten, 


Te = Tie YES: 


3.2.2 Hybrid pipe formulation 


To get a more accurate primal problem with still acceptable computational perfor- 
mance a combination of a pipe model with a static transport delay, calculated using 
the node method (Section 2.1.2), and a flexible transport delay, modeled using first- 
order differential equations (as the model presented in Section 3.2.1), is introduced in 
the following. This concept is based on [BHW17]. Because in real-world applications 
there is always an upper bound on mass flow and velocity, there is a minimum of the 


transport delay. This static minimum transport delay T"" can be calculated based on 
pipe length L and upper bound of velocity v4 or mass flow mU using 
2 
d 
L PIEZ 
TIRE u~ L 9 : (3.40) 
u 


If the upper bound of mass flow m“ is used for calculation, pipe diameter d and 


density of the transport medium p are required as additional parameters. 
The overall transport delay T; in a time slot t is split into this static minimum part min 


and a flexible delay 7; lex which is modeled using first-order differential equations. 
Hence, we arrive at a hybrid pipe model as shown in Figure 3.4 and the overall time 
delay 7; is calculated with 


m = T" | f, (3.41) 


i ex 
qin J 


0 L 


Figure 3.4: Scheme of a pipe with hybrid pipe model 


Modeling the remaining variable delay 7; = with a similar scheme as shown in Section 
3.2.1 the maximum mass flow rn needs to be subtracted from the mass flow m; in 
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the first-order differential equation in (3.39) resulting in 
2 
d 
px ($) i 1 1 ( av — Ti?) =i av — Ti?) 
n my mu i,t+1 thse i-1,t i,t 


kyndL (1 1 
Cpr 


T =r) (T7 = T) Vi € S VE E€ Sten. (3:42) 
Compared to the finite-volumes model presented in Section 3.2.1, the inflow temper- 
ature T°, , for the first volume element i = 1 is here replaced with the pipe inflow 
temperature T!” delayed by the constant minimum transport delay t’””. If heat losses 
are considered we arrive at a formulation as in (2.15) and get 


min 


m = Tb 4 (TR Tem)" Ye SiH I 3.43 
i-1,t 7 t— "in t,ı1 = l. ( Á ) 


The average temperature T;,; of the last volume element n remains the outflow tem- 
perature T?“' of the pipe in time slot t. The resulting model formulation for a primal 
problem using the hybrid pipe model is presented in Model 3.3. 
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Figure 3.5: Simulation results comparing finite-volumes and hybrid pipe models with an exact modeling 
of the transport delay for a pipe with 13.5 km length and 0.7 m diameter. 100 volume elements 


are used for the finite volumes model. 50 volume elements and a maximum flow of 3 m/s are 
used with the hybrid pipe model. 


This hybrid pipe model enables a more accurate representation of the pipeline dy- 
namics than the finite-volumes model presented in Section 3.2.1, as the smoothing of 
the delay only acts on of lex The temperature propagation including losses can be ex- 
actly calculated for the static minimum delay t””. This advantage of the hybrid pipe 
model is illustration in a simulation in Figure 3.5. However, depending on the combi- 
nation of discretization time, minimum delay time and pipe length the advantage of 
the hybrid pipe model to the finite-volumes model (Section 3.2.1) varies. 
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Model 3.3: Primal problem using the hybrid pipe model formulation 
Objective function 


Lb b . 
min ), (prices! YY — price Sell . pi) 


tes; 


+ D 2 ion nit + GR Dip + cian i ait) 


teS; i€Sg 


Generation and load constraints 


nit PE < pit < nig PH Vi € Sg, Yt € Sı 
nip QF < qis < nip QF Vi € Sg, YE € St 
b 
pimini- T Pit + p; uy el WEE S; 
IE Sg 
Fit = Cpi, (7; EN DR Vi € Sg U Sa, Yt € Sı 


The constraints for the feasible regions of the CHPs need to be added matching 
their type (e.g. (2.3) for CHPs with extraction condensing turbine). 


Node balances 
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3.3 Global optimization algorithm 


Given the dual problem presented in Section 3.1 and one of the primal problems pre- 
sented in Section 3.2 we can now combine them to a primal-dual global optimization 
algorithm for problems with variable-dependent time delays. For heating grids, there 
are non-convex bilinear terms besides the variable-dependent time delay. Therefore, 
to get a global optimization scheme for heating grids we extend the global optimiza- 
tion scheme of multiparametric disaggregation (see Section 2.2.2) with the dual and 
primal models presented in this chapter. 


Solve dual MILP problem 
Initialize primal non-linear problem 
with dual solution 


Solve primal problem (with fixed 
integers) using non-linear (local) 
solver 


decrease p and increase 
np (increase accuracy of 
dual problem) 


Check 
primal-dual 


gap 


Gap too big 


Figure 3.6: Iterative solution approach of multiparametric delay modeling 


This extension of the iterative scheme of multiparametric disaggregation is shown in 
Figure 3.6. The dual problem combines the outer approximations of multiparamet- 
ric disaggregation (Section 2.2.2) and multiparametric delay modeling (Section 3.1). 
Thus, an increase of the model precision requires a higher precision of approximation 
of the bilinear terms (by a decrease of p) and a higher precision of approximation 
of the variable-dependent delays (by an increase of np). One of the formulations pre- 
sented in Section 3.2 is used as primal problem modeling the variable-dependent time 
delay. 


The number of iterations needed to achieve an acceptable gap between the solutions 
of the primal and the dual problem depends on heating grid size and complexity. In 
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theory, with unlimited computational power and time the solution gap approaches 
zero. In reality, however, a gap will remain in most cases. For complex problems, a 
gap of about 1 % is usually accepted being one order of magnitude above the default 
gap of off-the-shelf MILP solvers (usually 0.1 %) [Gur16]. 


34 Summary of chapter 


This chapter presented “multiparametric delay modeling”, the first approach proving 
global optimality of its solution for optimization of district heating grids. “Multipara- 
metric delay modeling” is a deterministic global optimization scheme based on the 
decomposition of the original problem into a primal and a dual problem that are eas- 
ier to solve than the original problem. First, Section 3.1 introduced a novel dual model 
formulation to represent variable-dependent delays. It approximates the bounds of 
the pipe outflow temperature by introducing binary variables representing the trans- 
port time in every time slot. Thus, the resulting problem is a MILP problem which 
can be solved with a variety of efficient off-the-shelf solvers like [Gur16]. The out- 
flowing volume per time slot is split into sub-volumes to refine this approximation. 
Theorem 3.1 and Theorem 3.2 prove that this outer approximation becomes more and 
more accurate with increasing number of sub-volumes. 


In Section 3.2, two different primal models were introduced which can be solved with 
off-the-shelf NLP solvers using e.g. an interior point method [WB06]. In the first 
primal model, the pipe is split into a number of equal-sized finite volumes. For each 
volume a constant temperature and perfect mixing are assumed allowing to represent 
the thermal dynamics of the volume in one energy balance. The second primal model 
combines this finite volumes model with a static time delay leading to a hybrid pipe 
model. This hybrid pipe model allows a more accurate representation of the transport 
delay than the pure finite volumes approach of the first primal model. 


Finally, Section 3.3 established an extension of the global optimization algorithm of 
multiparametric disaggregation (Section 2.2.2), which enables to globally optimize the 
fully dynamic problem of heating grid optimization. For this, multiparametric disag- 
gregation for global optimization of bilinear terms is combined with “multiparamet- 
ric delay modeling” for global optimization of variable-dependent time delays. The 
dual problem presented in Model 3.1 and one of the primal problems (Model 3.2 or 
Model 3.3) are used in this global optimization scheme. 


4 Optimization with hybrid 
discrete-continuous time grid 


As the solution algorithm in the previous chapter is an iterative optimization scheme 
with multiple solver runs, solution times can be quite long. This chapter aims at 
finding an exact model which enables optimization in one solver run, promising faster 
solution times than an iterative solution approach. Recent advances in modeling in 
other areas help to achieve this goal. Here, we extend formulations from pipeline 
scheduling for refined, liquid petroleum products (as presented in Section 2.3) to 
become suitable for heating grid optimization. This approach was developed in a 
collaboration with Pedro Castro from University of Lisbon, Portugal and is published 
in [MC20]. It adapts the pipeline scheduling model presented in [Cas17]. 


As presented in Section 2.3, there are many similarities between pipeline scheduling 
for refinery products and optimization of heating grids considering the grid dynam- 
ics. Limiting the choice of supply temperature to a discrete set of temperature levels 
allows to use pipeline scheduling approaches for district heating grid optimization. In 
addition, a hybrid discrete-continuous time grid is introduced adapting the continu- 
ous time formulation of pipeline scheduling models to time discrete energy prices. 


4.1 Hybrid discrete-continuous time grid 


The choice of an appropriate time representation is crucial for scheduling models 
[HMB*14, MHI*15]. Looking at the literature review in Section 2.1.2, Section 2.1.4 
and Section 2.1.5, discrete time models are dominating past publications on heating 
grid optimization. However, for pipeline scheduling continuous time representations 
have shown to be very efficient to model flow dynamics (see Section 2.3). 


A combination of the two is presented in [WL18] where a hybrid discrete-continuous 
time grid models thermal dynamics of a district heating grid. In particular, fixed time 
points are used at every border between periods of constant demand and prices and 
combined with floating time points to account for the variability of the time delay. As 
the number of floating points between two fixed time points depends on the thermal 
grid setup, an upper bound of two per thermal storage tank in the system is proven in 
[WL18], but no model formulation suitable for optimization is presented. Therefore, 
the hybrid-discrete continuous time model for district heating grid optimization is 
introduced in the following as first presented in [MC20]. In this model thermal storage 
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tanks are not considered, but the dynamics of the district heating grid are used as 
thermal storage. Heuristically, one floating point is defined between two fixed time 
points based on optimal performance in a small case study. In systems with a larger 
number of pipes and nodes, a higher number of floating time points between two 
fixed time points might be needed, increasing the complexity of the formulation. 


Actual time 0:00 1:00 2:00 3:00 23:00 24:00 
Time grid ee ec = +H 
Fixed points 1 3 5 7 47 49 
Floating points 2 4 6 8 46 48 
are 
Variables ty litz t tst ty tg t46 t47 t4gt49 


Figure 4.1: Hybrid discrete-continuous time grid 


Figure 4.1 shows the hybrid discrete-continuous time grid. All odd-numbered time 


points are located at hour marks tfx, and are thus fixed time points t € sf “All 
even-numbered time points are floating time points allowed to move freely between 
the neighboring fixed time points. For one day or 24 hours we get a total of n; = 49 
time points (see Figure 4.1) and accordingly 48 time slots as every hour is split into 
two time slots by the floating time point. A big advantage of this time representation 
in comparison to a continuous time representation is that the electricity price price?! 
for each time slot is known a priori and a complicated assignment of variables to time 
slots as for example in [HH13] is not needed. The same applies to the heat demand 
qgemand which is forecasted with hourly changing values. 


Let us introduce non-negative real variables to model this hybrid discrete-continuous 
time grid. t; represents the actual time at time point t, whereas L; represents the 
duration of the time slot. 


t4ı = ti + Lt Vt € Sten (4.1) 


allows to link the two. All even time points are fixed time points tfx, which is enforced 
by 
het, esf. (4.2) 


All odd-numbered time points are floating time points t ¢ sf * which are allowed to 
vary between the neighboring time points in interval [t;_1, ¢;41]. This can be guaran- 
teed by limiting the time slot length L; to one hour, the distance between two fixed 
time points. As a time point is needed, whenever there is a new batch of water enter- 
ing or leaving a pipe and this does usually not happen at the hour marks, the floating 
time points allow an exact representation of varying in- and outflows. In some cases 
there will be no need for a floating time point between the hour marks e.g. when 
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there is no change in the in- or outflow of the pipe. This leads to dummy time slots 
with L;_ı = 0 or L; = 0 as the optimization chooses t+ equal to either t;—1 or tı 11. 


Remark: As explained in [MC20], choosing a continuous time grid where all hour 
marks need to be assigned to a time point and all time points are floating time points 
as proposed in [CHG09] does not lead to superior solutions. At least n; = 26 time 
points are needed to get a feasible problem, being one above the minimum number 
for one time slot per hour. But in comparison to the formulation using the hybrid 
discrete-continuous time grid, this continuous time formulation has a larger problem 
size, worse results and orders of magnitude longer solution times. 


4.2 Supply node 


Besides the hybrid discrete-continuous time grid another novelty in [MC20] is the 
discretization of temperature levels that enables to adapt the product-centric pipeline 
scheduling formulation [Cas17] to district heating grid optimization. The following 
formulations are introduced for one CHP, but adding more CHPs is straightforward 
by adding an additional index to all CHP variables and parameters. 


In comparison to the basic unit commitment model presented in Section 2.1.1 the bi- 
nary variable n; denoting the running state of CHP i is extended with one additional 
index l, being the temperature level currently produced by the CHP. Thus, for each 
CHP we now have a binary variable Xf, = 1 if in time slot t the CHP is producing 
water at temperature level /. As at most one temperature can be produced in a time 
slot, we get 

E XGI1 VEE Sten (4.3) 

leS; 


The set S; contains all possible temperature levels 1. When the CHP is producing water 
at temperature level I in time slot t, the volume Fo (in m?) enters the pipe, being a 
non-negative real variable with lower and upper bounds f! and f¥. In comparison 
to the lower and upper bounds Q} and QU on thermal energy generation of the 
CHP presented in Section 2.1.1, these are bounds on the volume injected by the CHP 
which can be calculated for example based on minimum and maximum mass flow or 
velocity: 

S XG SFG < fUXE, VIESWVIE Sten. (4.4) 


The change of inner energy of the transport medium at the producer (see (2.9)) needs 
to account for the different temperature levels ! as well. The thermal energy volume 
qt produced in time slot t is calculated based on supply temperature levels T; “pply 
and return temperature T’*"" (assumed to be constant) with 


1 
ge = BG (Teen) Vie Stene (45) 
1 
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As before, c, is the specific heat coefficient of the transport medium (water) and p is 
its density. 


The equations for the electric generation as well as the cost function remain the same 
as for the basic unit commitment model (Section 2.1.1), when both time slots within 
one hour are assigned to the parameters of this hour. The equations linking electric 
and thermal output of the CHP, however, need to consider the variable length of 
the time slots. Thus, formulation (2.3) for the feasible region of CHP operation with 
extraction condensing turbine is not suitable. To match the real feasible operation 
region (Figure 2.1) the inequalities (2.3) are adapted to 


pu Sau- Vt € St, 
a > az- Li — Body VEE St, (4.6) 
Su > —a3- Li Bon VEE S; 


with parameters «1, &2, &3, 61, 2 and 3 adjusting slope and offset of the region 
boundaries. 


4.3 Demand node 


Similar to the supply node, let us introduce a binary variable xP , being 1, if the water 
arriving at the demand node is at supply temperature level I in cae slot t. Due to the 
transport delay in the pipe, this is not necessarily the supply temperature level at the 
CHP in time slot t. At most one temperature level ! can arrive at the demand node in 


one time slot t: 
EXB S1 VEE Sten. (4.7) 
l 


The volume of water at supply temperature level ! consumed by the demand node FR 
(in m?) in time slot t is bound due to limitations in mass flow and velocity. Similar to 
the volume produced by the CHP we get 


SEXP, < FR < fX, VIESWVEE Sten. (4.8) 


As the mass flow to the demand node is limited to interval [m], m], the operating 
volume flow needs to stay always within these bounds as well. Without an explicit 
variable for the mass flow, we can ensure these bounds with 


1-)OX VEE Sten. (4.9) 
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This is possible, because mass and volume flow are connected via the density p of the 
transport medium. For this equation it is important that no time slot is longer than 1 
hour. 


The energy balance for the demand node cooling water from supply temperature 


level T;"?? 'Y to return temperature T’°“"" is explained by the change of inner energy 
of the transport medium based on (2.9). As for the CHP, all different possible tem- 
perature levels ! need to be considered. In addition, the thermal demand needs to be 
scaled according to time slot length L; to guarantee a continuous delivery of thermal 
energy: 


gaemand Re DER: P- Cp (7, poupply — return) VEE Sten: (4.10) 


If thermal losses are to be considered, the supply temperature level To PIY in the pre- 
vious equation needs to be adjusted to represent the temperature aioe along the pipe. 
As the supply temperature is the major influence on thermal losses (2.15) each supply 
temperature level at the demand node must be corrected individually by introducing 
new parameters representing the reduced supply temperatures per temperature level 


at the demand node. 


4.4 Pipe representation 


To model the flow dynamics in a pipe, we introduce the non-negative real variable V) ; 
which represents the volume of water at temperature level / inside the pipe in time 
slot t. It is assumed that there is only one volume of temperature level | inside the 
pipe in every time slot t. In Figure 4.2 all scenarios changing this volume are shown: 
If the CHP is feeding a volume ay at temperature level / into the pipe in time slot 
t, the volume V;; at temperature level / inside the pipe increases as shown in Figure 
4.2 a). If water at temperature level I arrives at the demand node, the volume V}; at 
temperature level / inside the pipe decreases by FR as shown in Figure 4.2 c). Figure 
4.2 b) shows a situation, where the volume V/, at temperature level I inside the pipe 
does not change in time slot t. This is possible as well if V;; = 0 or Vj; = v ina time 
slot. The respective volume balance for the volume of water at a certain temperature 
level I inside the pipe in time slot t is thus 


Vie = ofl + (Vier +F- Abs) lo: YESYES. (11) 


Here OP | pad is the initial volume of water at temperature level / inside the pipe. 


As water is non-compressible, the sum of all volumes inside a pipe needs to equal the 
overall pipe volume v in all time slots t: 


\V,=v WEE Sp. (4.12) 
l 
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Figure 4.2: Volume V;; at temperature level ! a) flowing in, b) being in and c) flowing out of a pipe with 
overall volume v. Note that coordinates LC); and RC); are volume coordinates, not distance 
coordinates. 


As shown in Figure 4.2, the non-negative real variables LC;; and RC); describe the 
left and right coordinates of volume V; having temperature level / in time slot tł. Note 
that these coordinates are not distance coordinates in m but volume coordinates in m? 
describing the position of water at temperature level | inside the pipe with respect to 
the overall pipe volume v. Accordingly, they are limited by pipe volume v: 


IC), <V VI E S, VEE St, (4.13) 
RC;; <v VI € 51, Vt E St. (4.14) 


The coordinates LC;; and RC); need to be directly left and right of the volume at 
temperature level / and, thus, have a difference of 


RCs — LC = Ve YLE Sy, VE € Sy. (4.15) 


To complete the model some additional variables need to be introduced. Variables 

LC and RCPS describe the coordinates of water at temperature level ! at the 
end of the last time slot t — 1. These coordinates are different from the coordinates 
at the beginning of the time slot if water at a different temperature level entered or 
left the pipe. To link the coordinates at the end of a time slot Loe and RG with 
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coordinates at the beginning of the time slot LC); and RC) +, we use 


Lc = LC + VES VI € S}, Vt € Sten, (4.16) 
Lt ß It i 
a 
RCHA = RC + YEP VI € S1, YE € Stem. Cp 
Lt , It : 
VA 


We propagate coordinates Loe and RO, at the end of the previous time slot t — 1 
to the coordinates LC); and RC); at the beginning of time slot t using 


LC} = Ie? |i + (taa 1- ZC) lt>1 VI E S),Vt € Si, (4.18) 


RCi = rellı=ı + (RCP — ZCrp) la WE SYES (4.19) 


In addition, these equations initialize the coordinates at the beginning of the optimiza- 
tion horizon to the initial values of the coordinates Ic? and rc? for every temperature 
level I. 


Variable ZC); is needed to reset the coordinates if water at temperature level I is not 
inside the pipe, such that it can enter the pipe. Hence, it is linked with the binary 
variable X, indicating if a temperature level I is currently inside the pipe: 


AEE (1 = Xi) VIE Sp YEE Sion. (4.20) 
Accordingly, the volume V;; of water at temperature level / inside the pipe in time 


slot t as well as the corresponding left and right coordinates LC); and RC;; need to 
be zero, if there is no water at temperature level | inside the pipe (XP, = 0). 


Vit < v- Xf, VIE S, Yt € Sısı (4.21) 

LCi; < v- XP VI € S Yt € Sısı (4.22) 
, Lt 

RC pose VIE Si, Yt € Sis (4.23) 
, Lt 


Finally, the conditions describing when water at a temperature level ! can enter or 
leave the pipe need to be defined. As it is assumed that only one volume of a tem- 
perature level is possible inside the pipe, the left coordinate LC), for this temperature 
level I has to equal zero if water at temperature level | should enter the pipe (Xf, =1) 
and we get 


LC <0: (1- XG) WES), VEE Sten (4.24) 


This equation covers both cases: if there is no water at temperature level ! inside the 
pipe (XP, = 0 and V;; = 0) as well as if water at temperature level / is inside the 


pipe (xP , = land 0 < V); < v) and was entering in the previous time slot as well 
(RC: = Vit). 
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For water at temperature level | leaving the pipe (XP, = 1) we need to ensure that its 


right coordinate RC); is at the most right place and, thus, equal to the pipe volume 
v: 
RC > v: xp, VI E S, VE E€ Stene (4.25) 


4.5 Overview of hybrid time grid model formulation 


Combining all equations presented in this chapter can replace the heat balance with- 
out thermal dynamics (2.5) in the unit commitment problem presented in Section 2.1.1. 
This results in a MILP formulation for the scheduling of CHPs considering heating 
grid dynamics. This resulting model can be solved directly with off-the-shelf MILP 
solvers and an overview of this model is given in Model 4.1. 


Model 4.1: Hybrid discrete-continuous time grid model 
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Feasible region of CHP with extraction condensing turbine: 
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Pipe representation 
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ZC < fu: (1- Xf) YI E Sp Yt € Sis 
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4.6 Summary of chapter 


The model presented in this chapter and published in [MC20] is based on latest find- 
ings in pipeline scheduling which are explained in Section 2.3. The supply temper- 
ature was discretized to adapt these scheduling models, that have been designed for 
scheduling transport of refined petroleum products via pipes, to heating grids. Thus, 
only a discrete set of supply temperatures is allowed. Of course, this discretization 
reduces the possible solution space possibly leading to worse results. However, if 
enough temperature levels are introduced (e.g. every 5 Kelvin), the solution space 
should remain big enough that this influence remains small. The original pipeline 
scheduling model uses a continuous time representation [Cas17]. In contrast, our 
model uses a hybrid discrete-continuous time grid (see Section 4.1) that enables an 
efficient integration of energy prices and load forecasts that usually vary at discrete 
time points (e.g. every hour or quarter hour). These two extensions, the discretization 
of supply temperature and the discrete-continuous time grid, result in an optimiza- 
tion model for heating grids which enables an accurate representation of transport 
times with floating time points, while staying solvable with off-the-shelf solvers as 
the resulting problem is a MILP and non-linearities are avoided. 


The hybrid discrete-continuous time grid is unique in modeling of heating grids, 
as approaches in literature usually use a discrete time representation. The pipeline 
scheduling models presented in [Cas17] allow to consider meshed grids and reversible 
flows. As the approach presented in this chapter is an extension of these models, the 
additions allowing meshed grids and reversible flows could be applied for heating 
grids as well. Consideration of pumping cost and an accurate representation of heat 
losses will be difficult, as using a volume-centric model formulation there is no explicit 
variable for mass flow or velocity. However, an integration of approximate heat losses 
depending on the current supply temperature level is possible by introducing one new 
parameter per temperature level which represents the reduced supply temperatures 
at the demand node. 


The hybrid time grid approach requires detailed knowledge of grid topology and pipe 
parameters (e.g. length, diameter, heat loss coefficient). Thus, it seems not suitable 
for distribution grids with many small pipes, but very suitable for heat transport 
networks with few large pipes. With the innovative problem formulation it allows an 
exact calculation of flow duration while resulting in a well to solve MILP problem. 
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One of the major drawbacks of the optimization approaches considering heating grid 
dynamics discussed in state of the art as well as in the previous chapters is that they 
all require a detailed grid topology. For real-world use of these models, this leads to a 
modeling error, because very often not all pipe parameters are known. Furthermore, 
this requires a high engineering effort, as the heating grid needs to be modeled in de- 
tail before an aggregation strategy or another simplification approach allows to build 
suitable optimization models. In addition, most methods use either an extensive itera- 
tive or sequential solution algorithm or strong assumptions (like constant mass flow). 


Grof uses a different approach by proposing a function f directly linking thermal 


energy stored in the heating grid qo with current and past supply temperatures of 


generators T; "" W and the current overall load factor geemand y gdemand max. [Gro12]: 


demand 


charge _ supply „supply „supply „supply q 
qt “(x ‚T_ı ‚12 /T 3 rt 


qdemand max 


VEE S; 61) 


Using such a function as in (5.1) enables building a very fast and efficient optimization 
model if added to a unit commitment and dispatch model without consideration of 
grid dynamics (Section 2.1.1). In [Gro12] this function is built using a linear regression 
trained with simulated data. Hence, the engineering effort needed to build a detailed 
heating grid simulation model remains. Using linear regression mainly shifts compu- 
tational effort from optimization runs to optimization model parameterization. 


In the following, a similar function as (5.1) is derived. However, to reduce engineering 
efforts for real-world applications, a more direct way to find this function is proposed. 
Section 5.1 recalls the storage effect of thermal dynamics. Section 5.2 introduces an 
outer approximation of the thermal energy stored in the heating grid which leads to 
the delay matrix approach in Section 5.3. This approach was developed in a master 
thesis [Hail8] and is published in [MHH19]. It allows to consider thermal dynamics 
of heating grids with all thermal generation in one place. The integration of thermal 
losses presented in Section 5.3.1 and a detailed real-world case study using this ap- 
proach are described in [MKV*19]. Section 5.3.2 presents an extension to multiple 
CHPs in different locations that was developed in an internship [Lor19]. 
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5.1 Dynamics of a heating grid with multiple loads 


To recall the main principles of thermal dynamics and the thermal storage effect in a 
heating grid we extend Example 2.2 with one producer and one load to Example 5.1 
with one producer and two loads following [MHH19]. 
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Figure 5.1: Storage behavior of a heating grid with two consumers. a) Supply temperatures at producer 
and loads. b) Heat production of producer and sum of heat consumption of consumers. 
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A situation with one producer and two loads is shown in Figure 5.1. Here the increase of 
supply temperature at the producer in time step 3 reaches the first consumer in time step 5 
and the second consumer in time step 7. Thus, the increase of the thermal energy output 
of the producer starts in time step 3 and is reduced to its original level in two steps. When 
the increased supply temperature reaches the first consumer in time step 5 this consumer 
reduces its mass flow. Thus, the overall mass flow is reduced, but the second consumer is 
still supplied with increased mass flow. This reduction is proportional to the mass flow or 
load share of consumer one. When the supply temperature increase reaches the second and 
last consumer, the thermal energy output of the producer is back to its original level and, 
hence, the charging of the grid ends. 

For a decrease of supply temperature at the producer in time step 9 we can observe an 
inversed behavior (similar as in Example 2.2). The thermal energy output of the producer 
is reduced immediately leading to a discharge of the grid. The thermal energy production 
gets back to its original level in two steps. When the first consumer is reached in time 
step 11, the producer increases its output with the now increased mass flow of consumer 
one. When the second consumer is reached in time step 13 the thermal energy output of the 
producer reaches its original level and discharging of the grid ends. Similar to Example 2.2 
thermal losses are represented in the reduced supply temperatures at the loads. 


5.2 Outer approximation of grid storage behavior 


To integrate the thermal storage behavior of a heating grid explained in the previ- 
ous chapter into a basic scheduling problem for CHPs (see Section 2.1.1), we try to 
extend the standard models for electric battery storage. Thus, we introduce one non- 


negative, real variable representing the state of charge SOCH*“! in time slot t and one 


real variable representing the thermal energy charged or discharged He in time 


slot t. Using these variables, the dynamic behavior of the state of charge is described 
with 
socteat — sochet + se re 5S. (5.2) 


To integrate thermal energy storage into the basic unit commitment problem (Section 
2.1.1), we adapt the thermal energy balance (2.5) to 


ge Vo gl tg er es, (5.3) 


iCSe i€S4 


For a battery, the state of charge SOC!“ is limited by the capacity limits of the storage 
and the energy of charge or discharge grs is limited by the maximum charging and 
maximum discharging power. For dedicated thermal storage like storage tanks these 
limits are usually known. However, it is more complicated to derive these limits for 


the inherent thermal storage of a heating grids. In the following we introduce such 
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limits suitable to approximate the thermal dynamics of a heating grid as we showed in 
[MHH19]. Looking at Example 2.2 and Example 5.1 there is a clear upper limit on the 
energy charged to as well as the energy discharged from the heating grid being linked 


to the maximum possible supply temperature difference (re —T, PP ) in time slot 


t. With the equation for the change of inner energy (2.9) we can relate this maximum 
temperature difference to a maximum energy volume charged to or discharged from 
the grid and get 

= cpt (TY — TPP) < GR < cpu (TUE) EES, 6A 
where TU is the maximum possible supply temperature, T upply is the originally 


planned (minimum) supply temperature in time slot t, 11; is the mass flow from the 
generation site in time slot t and c, is the specific heat capacity of water. 


The upper bound for the state of charge SOC}! of the thermal energy stored in the 
heating grid in time slot f is more complicated to derive. Assuming we know the 
maximum transport delay t’”* to the last consumer in the heating grid, we can sum 
up the maximum charging power over this time horizon to find an upper limit for the 
state of charge SOCH““ resulting in 


t 
0< sock < Y cpm, Ga = p VEE Sp. (5.5) 


tt=t— max 


However, with this upper bound of the state of charge SOCH“! we only get an outer 
approximation of the possible energy stored in the grid. In scenarios with multiple 
loads not all consumers are served with the maximum transport delay t'"“*. Accord- 
ingly, the thermal energy charged to the heating grid cannot reach the maximum level 
for the full duration t’*. In Example 5.1 this formulation overestimates the storage 
capabilities of the heating grid, as in time step 5 the supply temperature increase 
reaches the first consumer (consuming 50 % of the thermal load) and in the following 
charging the grid with thermal energy is reduced by half. With the representation 
above it is assumed that maximum charging is still possible. 


Thus, to derive a better representation of the thermal storage behavior of heating 
grids, the delay matrix approach is presented in the following section as published in 
[MHH19]. 


5.3 Delay matrix approach for co-located 
thermal generation 


From Figures 2.4 and 5.1 it seems to be possible to derive a function that directly links 


supply temperature at the producer T:"?? ly and thermal energy charged or discharged 
Pply P P t 87 8 8 


from the grid qe as in (5.1), if constant mass flow is assumed. 
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Let us investigate the temperature and thermal energy variations in Example 2.2. The 
temperature increase from time step 2 to time step 3 is 


suppl suppl 
qo arated (5.6) 


The mass flow rm does not change between these time steps, which leads to an in- 
creased output of thermal energy of the producer of 


l l 
epi (TM = Ta) (5.7) 


When the increased supply temperature reaches the consumer in time step 5 after 
a time delay T = 2, the thermal output of the producer is reduced by the same 


amount: ; j i i 
Supp!y SUPPLY: Supply SUP Ply 
epi TY TEES eee (Te ra (5.8) 


As for the outer-approximation of the storage behavior in Section 5.2 we introduce 


gae as thermal energy being charged to or discharged from the grid resulting in 


the energy balance 


D qian +. Ge fe qg” = = L qit VEE Si. (5.9) 


i€S4 i€Sg 


In Example 2.2 with a constant head load, changes in the heat output of the producer 


directly influence the thermal energy stored in the grid Ta Hence, with (5.7) and 
(5.8) we get 


charge _ charge ; supply supply 
a tep (7; -Ta 


t—T 


; l l 
— cpn (TPY — TY) VEE Sis, 610 
where T; "PP IY is the supply temperature at the producer in time step t, T is the trans- 
port time from producer to consumer and 11 is the constant mass flow. For time 
varying mass flow rn; considered to be a parameter, (5.10) becomes 


genarge = rats +Cp (rin T? upply _ tis E) 
= (= de TER re) Vt © Sis. (6.11) 


However, the time varying transport delay coming with time varying mass flow is 
not represented correctly in the above equation as there is only one transport delay T 
possible. To cope with this problem, we introduce a matrix M allowing flexible trans- 
port delays. Every row l in this matrix represents the start of a temperature increase in 
time step l. This allows to have different transport delays throughout the optimization 
horizon in line with the varying mass flows 1i;. Rows | in matrix M have non-zero 
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Table 5.1: Delay matrix M [I,t] for Example 5.1 


Na 28 ae 2 


0 
0 
0 
0 
0 


entries in columns t only if an increase of consumption in time step l reaches a con- 
sumer in time step t. The value of this non-zero entry is the share of consumption of 
the respective consumer. 


In Table 5.1 the matrix for Example 5.1 is shown. As the transport delay to the first 
consumer is two time steps, the increase of supply temperature at the producer in 
time step | = 3 would reach the first consumer in time step t = 5. The second con- 
sumer is reached after a transport delay of 4 time steps. Thus, the increase of supply 
temperature at the first consumer in time step l = 3 reaches the second consumer in 
time step t = 7. As both consumers are having equal load share, the non-zero entries 
of delay matrix M in row 3 are 0.5 in column 5 and 0.5 in column 7. For all other time 
steps l the rows are filled accordingly. If the transport delay is not an integer number 
of time steps, the consumption share can be split up between two or more columns to 
achieve a better representation of the real-world delay. 


Adjusting the calculation of qe in (5.11) to integrate delay matrix M we get 


charge charge . suppl 5 suppl 
rt Cp Kar ml A 


nt nt 
— cp (Em Lim TP Y Mt]. marp) Yt E S1. (5.12) 
l=1 l=] 


As an alternative to this recursive calculation of qo that we introduced in [MHH19], 

an explicit formulation was developed in an internship [Lor19]. The energy charged 

to or discharged from a pipe is the difference between the thermal energy flowing 
out 


into qi" and out of gq?" the pipe: 
gp = git — gi = WEE Sp. (5.13) 


As the thermal energy feeding the pipe is explained with the inner energy balance in 
(2.9), we get 
gi" = cprin (a = pen VEE Sp. (5.14) 
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The thermal energy flowing out of the pipe is the inflowing thermal energy delayed 
with transport time Tt: 


gt = Gig = cyte (TT) VEE Se ae 


Hence, the energy charged into a pipe qe can be calculated with 


qE = cy (mi TPPP — tug TEY) Vie. (5.16) 
Thus, using delay matrix M[I, t] a similar formulation as in (5.12) is possible: 


Nt 
Pii Si (m F nee = L M |L t] - my - ae) VEE Sp. (5.17) 
l=1 


(5.12) or (5.17) in combination with (5.9) and delay matrix M (see Table 5.1 for Example 


supply 


5.1) give a direct link between supply temperature at the producer T; and thermal 


energy charged to or discharged from the heating grid gore. Thus, they can extend 
the optimization problem for production scheduling without consideration of grid 
dynamics in Section 2.1.1 to an optimization problem considering the storage effect 
induced by supply temperature changes. The consumption share of and transport 
delay to consumers as well as the mass flow at the producer must be known for this 
formulation. Model 5.1 shows an overview of this optimization problem at the end of 
the next section as it includes thermal losses introduced in this section, too. 


Remark: As many models in literature, this approach assumes fixed, given transport 
delays. It does not account for the change in transport delay when a consumer reduces 
its mass flow with an increased supply temperature arriving. Hence, there is a model 
error being larger for longer and higher temperature increases. This model error 
will reduce if variable mass flows are introduced. The approach presented in Section 
5.3.2 to integrate a second CHP can be a starting point for this extension as here two 
different mass flows are considered. 


5.3.1 Integration of thermal losses 


To consider the additional heat losses caused by the increase of supply temperature, 
we introduced the following approach in [MKV*19]. In real-world heating grids, 
there are usually only a few thermal generation points which are measured with high 
sampling rate and hundreds of loads normally not measured with this accuracy. Thus, 
thermal demand forecasts with high time resolution are usually predicted based on 
measurements of past generation at different ambient conditions. Accordingly, those 
load forecasts gee will include the thermal losses of the heating grid with the 
typical supply temperature control. However, increasing the supply temperature to 
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store thermal energy within the heating grid increases thermal losses. This should be 
considered as gee in (5.9). Recall (2.13) describing the heat losses of a pipe filled with 
a liquid at one temperature: 


qe — kyrdL (ete = re) . 
As ambient temperature T™? and pipe parameters do not change with an increase of 
temperature, the increased losses are proportional to the temperature inside the pipe 
Tinside_ Thus, we can introduce a loss factor dos, and using a time discrete formulation 


we get De 
qu = Kloss ' De (5.18) 


Now the question remains, how to best replace T/":ide with a formulation using the 


supply temperature at the producer T;"?? W. Heat losses first lead to a reduction of 
temperature of the transport medium. This does not effect the thermal energy output 
of the generation until this lowered temperature arrives at the consumer and the 
consumer increases its consumed mass flow. Using delay matrix M we can describe 
this delayed arrival by 


N 
qos = Nipgs * y M [L t] TpuPPly VEC St. (5.19) 
l=1 


In comparison to the dynamic thermal loss calculation in (2.15) this equation nei- 
ther depends on mass flow m+ nor on transport time 7. If T; “PP!Y is the pipe inflow 
temperature, it overestimates the thermal losses as in reality the temperature drops 
along the pipe. However, for typical operation scenarios this linearized loss calcula- 
tion should represent the thermal losses sufficiently well, as transport time Tt has a 
minor influence on thermal losses in comparison to inflow and ambient temperature. 
The influence of transport time 7; in (2.15) grows only for very slow velocity and thus 
high transport times. Figure 5.2 shows the thermal losses for different pipe diameters 
and different supply temperatures at different typical velocities in heating grids as 
calculated from (2.15). It illustrates the behavior of thermal losses described above, as 
the velocity is proportional to the mass flow and inverse to the transport time. 


Model 5.1 summarizes the resulting model formulation for the delay matrix approach 
considering additional thermal losses induced by the increased supply temperature. 
These additional thermal losses should be considered in the operations planning prob- 
lem, especially if a distribution heating grid with a high number of pipes with small 
diameters is optimized. For transmission grids, which have only pipes with large 
diameters, thermal losses are often insignificant. 
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Figure 5.2: Thermal losses in a pipe with 5 km length and a) 0.1 m or b) 0.7 m diameter depending on 


flow velocity for different inflow temperature levels in 10 K steps from 70 °C (lowest) to 130 °C 
(highest) 


Model 5.1: Delay matrix model for co-located CHPs with thermal losses 
Objective function 
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The constraints for the feasible regions of the CHPs need to be added matching 
their type (e.g. (2.3) for CHPs with extraction condensing turbine). 


5.3.2 Integration of distributed thermal generation 


To reduce thermal losses, heating grids are usually supplied at minimum possible 
temperature. As there is a limitation on velocity and flow in a pipe, this minimum 
temperature can be calculated based on (2.9). To enable integration of distributed 
thermal generation into the delay matrix approach assume in the following that all 
CHPs will always feed at minimum temperature and, thus, maximum flow. 


: 
pipe Ei I 
CHP1 load area 


Figure 5.3: A heating grid with a second smaller CHP 


To introduce additional CHPs, first, we look at the situation shown in Figure 5.3. 
Here, one long pipe is connecting one large CHP1 with a load area. A second smaller 
CHP2 is directly located at the load area feeding the load area without important 
transport delay. Assuming we always supply with minimum supply temperature and 
thus maximum flow, we get two possible flows in the pipe in every time slot: One 
higher flow if CHP2 is turned off and one slightly lower flow if CHP2 is running. 


If CHP2 is switched on (nz; = 1), it provides its maximum electric py and heat OF 
output: 


0 n24 = 0 

Pzt = í VEE St, (5.20) 
te n= 1 
0 n2t = 0 

12t = i VEE Si. (5.21) 
P n24=1 


At times when the thermal energy storage of the grid is not used (times of constant 
supply temperature), the thermal output of CHP1 q,; is reduced accordingly, leading 
to 


demand loss 
_ IN +4 n,=0 
qt = no _ „loss not = 1 VEE St. (5.22) 


Cht — 42, 
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Hence, the mass flow m; in the pipes from CHP1 to the load area varies with the 
commitment status nz; of CHP2: 


1 _ı „loss = 
cp: (Te"PPY return ) (a! N ) n2t = 0 


iy = VEE Sp. (5.23) 


1 demand , „loss _ == 
cp: (TE Ten (a; qt 91) N21 1 


Thus, there is a limited amount of possibilities of temperature propagation in the pipe 
depending on current and past commitment states of CHP2. Assuming a scenario 
with a time delay limited to variations between two and three time slots, all possible 
temperature propagations are shown in Figure 5.4. The running status of CHP2 nz, 
in past time slots is encoded with Comb; as shown in Table 5.2. 


Table 5.2: Running status of CHP2 in past time steps depending on Comb; 


i Comb; mot Naı M21-2 


1 000 0 0 0 
2 001 0 0 1 
3 010 0 1 0 
4 011 0 1 1 
5 100 1 0 0 
6 101 1 0 1 
7 110 1 1 0 
8 111 1 1 1 


For the situation shown in Figure 5.4, matrix M [I,t] has only two entries per row and 
(5.17) becomes 

giS = ey (sing TE TE Tr) VE ES. 6.24) 
With lengths d! and d? defined as illustrated in the top pipe in Figure 5.4, parameters 


sl and s? can be calculated with the propagation Ax; of the transport medium in a 
time slot t as follows: 


1_ A 

ee (5.25) 
d2 

2. 

er (5.26) 


There is one s! and one s? for every combination i of current and past commitment 
states of CHP2 and thus current and past mass flows. As on the right side of Figure 
5.4, the values of s! and s? can be calculated for every possible combination i of 
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000 0.6 0.4 


001 0.6 0.48 


010 0.52 0.57 


011 0.52 0.68 


100 0.43 0.4 


101 0.43 0.48 


110 0.32 0.56 


111 0.32 0.68 


Figure 5.4: Possibilities of temperature propagation in a pipe with a major CHP1 at the beginning and a 
smaller CHP2 at the end (close to the load area) for different running states of CHP2 (bright: 
off, dark: on). The blue arrow indicates the flow direction. 
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current and past commitment decisions for CHP2 (Comb;). Introducing y+; as a binary 
variable being one if and only if case i is active, (5.24) becomes 


charge _ 


supply 
qi = cph: T, 


= Cp Luh orna + TEPPI 82. tya. 22) YEE Sh. (5.27) 


The mass flow m; in the pipe in time slot t depends on the combination of current 


and past commitment decisions of CHP2, too. As for every case i we know this 


1 
RER of RR and past mass flows, we can introduce new parameters s” 


and s” i ? combining s! ; and s2 ; with the corresponding mass flow: 
mi eh ite VIE Se (5.28) 


gt = s2. tapo VIE Se (5.29) 


Sc denotes the set of all possible combinations to run or not run the CHPs in current 
and past time slots. 


With these parameters, (5.27) becomes 


i 


genarge =Cy (m ; ae Ln ” (s mA en a ze 2) 
Vt €S;. (5.30) 


This equation still contains a multiplication of one binary variable (y; ;) with one real 


variable (TERP Y or To yy, Several MILP solvers, for example Gurobi [Gur16], are 
able to handle such equations. If the chosen solver does not support this type of 
equation, bigM constraints or generalized disjunctive programming can be used to 
formulate this equation as MILP [RG94, GBSt 12]. 


Of course, for every time slot t only one case i representing a combination of current 
and past mass flows can be selected. 


8 
\yi=1 WES: (5.31) 


The following equations ensure that y;; = 1 only if the correct combination of current 
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and past mass flows is active. Matching Table 5.2 we get 


nat +H not1+not-2t¥r2l Vie St 

not + nzt1 + (1— n22) +4221 VEE St 

nat + (1— not-1) + not-2 +yı3 2 1 Vt E St 

nat + (1— nzt-1) + (1 — n2t-2) tyta >21 VEES: 


(5.32) 


(=n) + Ss (1 nz2) +y >21 Mes. 


With Comb; as shown in the table in Figure 5.4, (5.32) can be reformulated more 
generally as 


3 
So nz 4-41 (1 — Comb; u) + (1 — not—t41) «Combi + yz; > 1, 
tt=1 


Vt € St,Vi € So. (5.33) 


(5.30), (5.31) and (5.33) allow to expand the delay matrix approach with a second CHP 
if they replace (5.12) or (5.17). Model 5.2 presents the resulting MILP model. If the 


mass flow varies more importantly with the switching of CHP2 as in the example 


shown in Figure 5.4, additional parameters sl, s?, 5 ... as well as more cases i are 


needed. (5.30), (5.31) and (5.33) need to be adapted to match this additional number 
of cases and parameters. Similar changes are needed if more CHPs are added to the 
optimization. However, the general concept remains: 


e enumerate all possible running state combinations of CHPs (Comb;) 
e for each state calculate parameters sl, s?, s$, u 


e adapt (5.30), (5.31) and (5.33) to the increased number of combinations 


Model 5.2: Delay matrix model for two CHPs considering thermal losses 


Objective function 
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Equality Constraints 
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The constraints for the feasible regions of the CHPs need to be added matching 
their type (e.g. (2.3) for CHPs with extraction condensing turbine). 


5.4 Summary of chapter 


This chapter presented the delay matrix approach, which offers a compact MILP for- 
mulation and, thus, promises fast optimization runs. Section 5.1 recalled the thermal 
dynamics of a heating grid and Section 5.2 showed how these dynamics are related 
with batteries and thermal storage tanks. Based on these findings, Section 5.3 pre- 
sented the delay matrix approach that summarizes the dynamics of the heating grid 
in a delay matrix M which only requires share of consumption of and transport time 
to loads as input parameters. Thus, it can easily be parameterized using measure- 
ment data which allows applications to real-world heating grids as we showed in 
[MKV119]. This is a major benefit in comparison to other optimization models as 
they usually require data of all pipes in the heating grid for parameterization (see 
Section 2.1.4 and Section 2.1.5). Section 5.3.1 presented an extension of the delay ma- 
trix approach that considers the additional thermal losses induced by storing thermal 
energy in the heating grid. As the thermal losses mainly depend on the transport 
medium temperature, a linear loss factor was introduced which can be estimated for 
real-world grids using the grid topology as in [MKV*‘19]. As many approaches in 
the literature, the delay matrix approach assumes constant transport delays. How- 
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ever, multiple transport delays can be considered as shown in Section 5.3.2. Here, the 
delay matrix approach for CHPs at one site was extended to a case with one major 
CHP storing thermal energy in the heating grid and one minor CHP at another loca- 
tion which is not storing thermal energy in the heating grid. To integrate this second 
CHP, different cases for all possible mass flow combinations are introduced in the for- 
mulation. This concept can be used to allow a flexible transport delay in the original 
delay matrix approach or it can be extended to represent more generation sites. 


6 Case studies 


In this chapter, the modeling and optimization approaches for district heating grids 
developed in Chapter 3, Chapter 4 and Chapter 5 are applied to different scenarios. 
In Section 6.1 two small district heating grid cases are introduced and used to com- 
pare the proposed optimization approaches as well as one approach from literature 
in terms of performance and solution quality. The comparison of the global optimiza- 
tion approach to the sequential approach from literature was published in [MLH19] 
including an evaluation of the different primal problems. In Section 6.2 the delay ma- 
trix optimization approach of Chapter 5 is applied to the district heating grid of the 
city of Kiel, Germany. This detailed real-world case study shows potential benefits for 
different generation setups and is published in [MKV* 19]. 


6.1 Small heating grid examples 


This section introduces to two small example heating grids, one with one CHP and 
one with two CHPs, which follow [MLH19]. For these two example grids, we present 
results for all optimization approaches introduced in this thesis and compare them 
with results of a sequential optimization proposed in literature. A detailed discus- 
sion of similarities and differences in performance of the algorithms concludes this 
section. 


6.1.1 Example cases 


Both example cases use the same electricity price profile and the same thermal and 
electric energy consumption profiles. For the electricity market it is assumed that 
electric energy can be sold and purchased from the day-ahead electricity market with 
a small premium for purchasing. As winter days are more promising for storing 
thermal energy in a heating grid, the electricity price profile is taken from EPEX SPOT 
DE/AT [EPE] for November 15th, 2017, a winter day with typical price variations, and 
shown in Figure 6.1. 


The thermal load profile in both scenarios is derived from a typical thermal demand 
for a November day of a German city at ambient air temperatures of 5 °C to 10 °C. It 
is shown in Figure 6.2. For the electric demand profile a scaled standard load profile 
for households shown in Figure 6.3 is used. 
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Figure 6.1: Electricity price profile (EPEX SPOT DE/AT, November 15th, 2017). Source: [EPE] 
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Figure 6.2: Thermal load profile for a November day with ambient air temperatures of 5 °C to 10 °C 
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Figure 6.3: Electric load profile (scaled standard load profile) 
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One CHP case 


Dee, 2 = 
pipe: 10 km length EB I 


CHP1 0.7mdiameter load area 


Figure 6.4: Example case with one CHP and one load area 


The first scenario considers the most simple heating grid possible: one load area fed by 
one combined heat-and-power generation plant (CHP1) via a pipe. Pipe dimensions 
are inspired by the district heating pipeline connecting Mannheim and Heidelberg 
[KC13] with a length of 10 km (length of Mannheim to Heidelberg connection is 
13.5 km) and a diameter of 0.7 m as shown in Figure 6.4. CHP1 is assumed to have 
a maximum output of 400 MW thermal and 500 MW electric power being generated 
with an extraction condensing turbine. Thus, the power-to-heat ratio is variable while 
being limited by the feasible operation region shown in Figure 2.1. CHP1 can heat 
water to a maximum supply temperature of 130 °C, but lower supply temperatures 
are possible as well. With (2.6) a linear cost function for operating CHP1 is assumed. 


Two CHPs case 


N CHP2 
pipe: 10 km length = 
CHP1 0.7 m diameter load area 


Figure 6.5: Example case with two CHPs and one load area 


The second scenario adds a second CHP to the first scenario. This CHP2 is located 
directly at the load area as shown in Figure 6.5 feeding the load area without a trans- 
port delay. CHP2 has typical characteristics of a gas turbine and can produce up to 
20 MW of thermal and 20 MW of electric power with a fixed power-to-heat ratio (be- 
ing 1). It heats water to a constant temperature of 110 °C. This setup is chosen as 
it is quite common to have one major heat supplier outside the main consumption 
area and additional smaller heat sources closer to the load area. Gas boilers or gas 
turbines are common choices for this additional heat sources due to comparably low 
investment cost, fast start-up and compact size. The characteristics of CHP1, the pipe 
and the load area remain the same as in the first scenario. A linear cost function is 


4 Heating grid pipes are under high pressure. Thus, water inside the pipes is liquid at 130 °C. 
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assumed for both CHPs. Table 6.1 shows a summary of the parameter differences of 
the CHPs. 


Table 6.1: Main parameters of the CHPs in the example cases 


CHP1 CHP2 
Electric output in MW < 500 <20 
Thermal output in MW < 400 <20 
Supply temperature in°C < 130 = 110 
Power to heat ratio see Figure 2.1 constant 


6.1.2 Global optimization results with 
multiparametric delay modeling 


In the following, the one CHP and the two CHPs cases presented in the previous sec- 
tion are optimized to global optimality. Multiparametric disaggregation (Section 2.2.2) 
is used to model bilinear terms and multiparametric delay modeling (Chapter 3) is 
used to model the variable-dependent time delay in the dual problems leading to the 
MILP formulation shown in Model 3.1. The primal problems are solved with a local 
NLP solver with direct implementation of the bilinear terms and using the formula- 
tions presented in Section 3.2 to model the variable-dependent time delays. For the 
primal problem using the finite-volumes model (Model 3.2), 500 equally sized volume 
elements are used to describe the pipe. Using the primal problem with the hybrid pipe 
model (Model 3.3), again 500 elements for the remaining flexibility are used. Thermal 
losses are considered in both the dual as well as the primal problems. 


As computational effort for the global optimization approaches is high, all global 
optimizations were run on an Intel Xeon CPU E5-2650 v3 with 2.3 GHz, 16 cores 
and 16 GB RAM using the modeling environment Julia/JuMP [Ran16], MILP solver 
Gurobi 8.1 [Gur16] for the dual problems and NLP solver IPOPT [WB06] to locally 
solve the primal problems. This case study using the global optimization algorithm is 
published as part of [MLH19]. 


One CHP case 


The result of the global optimization of the one CHP case using the global optimiza- 
tion scheme with the finite-volumes model as primal problem is shown in Figure 6.6. 
This result is the solution of the primal problem after 4 iterations having an overall 
objective value of 291,235.7€. The lower bound of the dual problem is 288,371.3 € 
which also provides a lower bound of the original problem. Hence, there is a global 
optimality gap of 1 %. The result has been achieved in 146.19 s. 
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Figure 6.6: Results of global optimization for 1 CHP scenario with finite-volumes model as primal problem. 


a) Supply temperature at CHP and at load. b) Thermal output of CHP in comparison to thermal 
load. c) Electric output of CHP, electric load and interactions with electricity market. 
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Looking at the electric production in Figure 6.6 c), it can be noted that in the early 
morning hours (hour 0 to hour 6) the electric output of the CHP varies a lot. Here, the 
price of electricity is approximately at the marginal cost of production. Thus, in one 
hour it is profitable to produce electricity whereas in the next hour it is not profitable 
to produce electricity. Running such big ramps in output of a coal fired power plant 
is usually not wanted. Therefore, in real-world operations one might choose either 
producing as much as possible or producing as little as possible in this time. Starting 
from hour 6 the CHP produces as much electricity as possible, as electricity prices are 
above the marginal cost of production. 


In Figure 6.6 a) and Figure 6.6 b) it can be noted that there are major drops in thermal 
output and supply temperature at hours 8, 12 and 18. From the electricity price 
profile (Figure 6.1) we note that these are times of higher electricity prices compared 
to the previous hours. Hence, here the supply temperature at the CHP is lowered to 
reduce the thermal output of the CHP enabling a higher electric output within the 
feasible region of the CHP shown in Figure 2.1. As reducing the supply temperature 
discharges the thermal storage capacity of the heating grid, there needs to be an 
increase of supply temperature beforehand to store energy within the heating grid. 
This increase is usually done as late as possible to decrease thermal losses. This can 
be seen in hours 10 and 11 directly before high-price hour 12, where an important 
increase of supply temperature charges thermal energy into the heating grid. At the 
end of the time horizon in hours 21, 23 and 24 after the last price peaks, the supply 
temperature remains at the lowest level, discharging the thermal energy stored in 
the heating grid. This is reasonable as the model formulation does not account for 
thermal energy stored in the heating grid at the end of the time horizon. 


However, an unexpected issue can be observed in Figure 6.6 a): The supply temper- 
ature is reduced to its minimum level at 60 °C several times within the optimization 
horizon. This is not feasible in real-world, as mass flow demand would exceed the 
maximum possible mass flow in a pipe if this temperature arrived at the load. How- 
ever, with the finite-volumes pipe model these low temperatures at the beginning of 
the pipe never reach the consumer due to the smoothing effect of the finite-volumes 
model. Thus, as shown in Figure 6.6 a) the temperature at the load stays within ac- 
ceptable bounds. Due to the same reason the temperature at the load does not look 
like a delayed version of the temperature at the CHP in Figure 6.6 a), which it is in 
real-world. 


The hybrid pipe model presented in Section 3.2.2 should allow for a more accurate 
representation of the temperature propagation. In Figure 6.7, the solution of the global 
optimization scheme with the hybrid pipe model as primal problem is shown. This 
result is the feasible solution of the primal problem after 4 iterations having an overall 
cost of 294,386.0 €. The corresponding global lower bound is 288,265.1 € leading to a 
global optimality gap of 2.1 %. This result was achieved in 125 s. 


In Figure 6.7 a) it can be seen that the supply temperature profile at the load follows 
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Figure 6.7: Results of global optimization for 1 CHP scenario with hybrid pipe model as primal problem. 
a) Supply temperature at CHP and at load. b) Thermal output of CHP in comparison to thermal 
load. c) Electric output of CHP, electric load and interactions with electricity market. 
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the supply temperature profile at the CHP a lot better than for the finite-volumes 
pipe model (Figure 6.6 a)). The electric output of the CHP shown in Figure 6.7 c) 
shows a similar behavior as before with a highly variable profile in the early morning 
hours up to hour 6 due to the electricity prices swinging around the variable cost 
of production. After hour 6, the electric output is at its maximum. As shown in 
Figure 6.7 b) the thermal storage effects are less important than with the previous 
model. As before, there is some thermal discharge of the heating grid (CHP output 
below load) around hours 8 and 18, and some charging of the heating grid with 
thermal energy in the hours before (CHP output above load). However, the amount 
of thermal energy charged to and discharged from the grid is less important than with 
the finite-volumes model. This can be explained with the reduced amplitude of the 
temperature changes at the CHP as smaller temperature changes directly lead to less 
thermal energy stored in the heating grid. 


Two CHPs case 


Figure 6.8 shows the result of the global optimization scheme (Chapter 3) using the 
finite-volumes model as primal problem (Section 3.2.1) for the two CHPs case pre- 
sented in Section 6.1.1. This result was achieved in 4 iterations and 1,637.0 s. It has an 
objective value of 348,252.9 € and a lower bound of 330,871.1 €. The global optimality 
gap is at 4.99 % accordingly. The solution time for solving the dual MILP was limited 
to 600 s in every iteration to limit the time spent solving one iteration. 


The generation of electric energy using CHP1 turns profitable in hour 7. Accordingly, 
the electric output of CHP1 is increased from minimum to maximum possible level 
as shown in Figure 6.8 c). CHP2 is switched on in hour 3 and from hours 7 to 21, 
being the higher price hours. This is reasonable as in contrast to CHP1 for CHP2 a 
higher thermal output means a higher electric output. Thermal energy is stored in the 
heating grid (generation of thermal energy above consumption) in the morning hours 
until hour 6, in hour 11, in hours 13 to 15 and in hour 22 as shown in Figure 6.8 b). In 
hours 8 to 10, hour 12, hours 18 to 19, hour 21 and starting hour 23 thermal energy is 
discharged from the heating grid (production of thermal energy below thermal load) 
with a reduction of supply temperature of CHP1 (see Figure 6.8 a)). As expected, 
these are hours with local electricity price peaks. As for the one CHP case, due to 
the smoothing effect of the finite-volumes model, the supply temperature at the load 
does not follow the supply temperature at CHP1 and the supply temperature at CHP1 
reduces to levels infeasible in real-world. 


The result of the global optimization of the two CHPs case with the more accurate hy- 
brid pipe model (Section 3.2.2) as primal problem is shown in Figure 6.9. Again with 
a limitation of the solution time of the dual MILPs to 600 s per iteration, this results 
was achieved in 6 iterations and 2,907.1 s. It has an objective value of 351,112.0€, a 
lower bound of 334,785.9 € and a gap of 4.6 %. 
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Figure 6.8: Results of global optimization for 2 CHPs scenario with finite-volumes model as primal prob- 
lem. a) Supply temperature at CHPs and at load. b) Thermal output of CHPs in comparison to 
thermal load. c) Electric output of CHPs, electric load and interactions with electricity market. 
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Figure 6.9: Results of global optimization for 2 CHPs scenario with hybrid pipe model as primal problem. 
a) Supply temperature at CHPs and at load. b) Thermal output of CHPs in comparison to 
thermal load. c) Electric output of CHPs, electric load and interactions with electricity market. 
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It can be seen that drops in supply temperature to unreasonably low levels occur less 
often with the hybrid pipe model as primal problem (Figure 6.9 a)) than with the 
finite-volumes model (Figure 6.8 a)). However, two large drops remain: one in hour 2 
and one in hour 11. The drop to 60 °C in hour 24 is reasonable for the optimization, as 
this temperature level does not arrive at the load within the optimization horizon. Ac- 
cordingly, thermal energy is discharged from the heating grid in hour 2 and hour 11. 
Additional discharging caused by less important and real-world feasible temperature 
drops occurs in hour 8, hour 13 and starting hour 17 as shown in Figure 6.9 b). Hence, 
the heating grid is discharged in higher price hour 8 and around high price hour 18. 
The discharge in hour 11, where the electricity price is lower than the surrounding 
time slots, seems to be counterintuitive, but might be induced by the starting values 
of the primal problem, which is only solved to local optimality, as the solver of the 
dual problem was stopped due to its time limitation. CHP2 runs in hour 3, hour 6, 
hours 8 to 12 and in hours 14 to 21 being the higher price times. Given the direct pos- 
itive link between thermal and electric output this is a reasonable choice. As shown 
in Figure 6.9 c) electric generation of CHP1 turns profitable in hour 7, as the electric 
output of CHP1 is increased to its maximum at this point in time. 


6.1.3 Optimization results with hybrid 
discrete-continuous time model 


The results for the optimization of the two scenarios using the hybrid discrete-con- 
tinuous time model of Chapter 4 are presented in the following. Heat losses are 
neglected, as for the studied cases with pipes of diameter 0.7 m they are of minor 
importance. 


These results were achieved on an Intel Silver 4110 CPU with 2.10 GHz and 16 GB 
RAM using the modeling environment Julia/JuMP [Ran16] and MILP solver Gurobi 
[Gur16]. 


One CHP case 


Figure 6.10 shows the optimization result using the hybrid time grid model for the 
one CHP case. To enable a direct comparison with the other optimization methods 
this result is averaged on hourly bases in Figure 6.11. Allowing eight temperature 
levels (90 °C, 95 °C, 100 °C, 105 °C, 110 °C, 115 °C, 120 °C and 130 °C) the solution 
with a MILP gap of 0.0513 % was achieved in 800 s resulting in an overall cost of 
292,565.9 €. 


In Figure 6.10 c) and Figure 6.11 c) we can observe that the electric output of the CHP 
in the morning hours from hour 0 to hour 6 shows the same volatility as with the 
global optimization in Section 6.1.2. This is caused by the same reason: The marginal 
cost of production is above the electricity price in one hour whereas in the next it 
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Figure 6.10: Results with exact timing for 1 CHP scenario with hybrid discrete-continuous time model. 


a) Supply temperature at CHP. b) Thermal output of CHP in comparison to thermal load. 
c) Electric output of CHP, electric load and interactions with electricity market. 
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Figure 6.11: Hourly averaged results for 1 CHP scenario with hybrid discrete-continuous time model. 
a) Supply temperature at CHP. b) Thermal output of CHP in comparison to thermal load. 
c) Electric output of CHP, electric load and interactions with electricity market. 
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is below. The CHP runs with maximum possible electric output after hour 6 where 
generation of electricity becomes profitable for the CHP. Looking at Figures 6.10 a), 
6.10 b), 6.11 a) and 6.11 b) the CHP runs with low supply temperatures and thus re- 
duced thermal output around hour 8 and hour 18. As explained in Section 6.1.2, this 
allows to increase the electric output during these hours of high electricity price. Ac- 
cordingly, the heating grid is charged with thermal energy by the supply temperature 
increases in hours 1 to 3 and hour 11 where electricity prices are lower than dur- 
ing the discharging times. This price difference allows to create economic benefits. 
As thermal losses are not accounted for, the time between charge and discharge of 
the heating grid is comparably long. With consideration of losses, this duration will 
reduce as every time slot with high supply temperature leads to increased thermal 
losses. 


To evaluate the accuracy of the hybrid time grid model, we compare the optimization 
result q1, of the thermal energy output of the CHP with the result of a simulation 
Gi, of the studied case. This simulation uses the node method (Section 2.1.2) with 
the optimization result of the supply temperature of the CHP and the thermal load as 
input. The comparison shows that the thermal generation of the CHP in the optimiza- 
tion result deviates with a root mean square deviation? of 29.431 from the simulation 
result. 


Two CHPs case 


To integrate CHP2 to the hybrid time grid approach, an energy balance between ther- 
mal demand, thermal output of the pipe and thermal generation of CHP2 is added. 
As before, eight temperature levels are allowed (90 °C, 95 °C, 100 °C, 105 °C, 110 °C, 
115 °C, 120 °C and 130 °C). Figure 6.12 shows the optimization result for the scenario 
with two CHPs being achieved in 800 s with an MILP gap of 0.0554 %. The objective 
value is 346,490 €. Figure 6.13 shows the same results but with hourly averaged values 
to enable a direct comparison with the results of the other optimization approaches. 


In the morning until hour 6 as well as from hour 11 to hour 17 the supply temperature 
of CHP1 is increased to store thermal energy in the pipe. Around hour 8 and hours 18 
to 22 the supply temperature is set to its minimum possible level, discharging thermal 
energy from the heating grid. Thus, in these hours the thermal output of CHP1 
is reduced and its electric output is increased to sell more or buy less electricity at 
higher price times. In the morning at hour 7, the operation of both CHPs turns 
profitable. Electric output of CHP1 is increased to its maximum possible level and 
CHP2 is switched on. In the evening, at hour 23 CHP2 is switched off as electric 
demand and electricity prices fall. 


5 The root mean square deviation is calculated with v Eti (gu - du) /n. 
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Figure 6.12: Results with exact timing for 2 CHPs scenario with hybrid discrete-continuous time model. 
a) Supply temperature at CHP1. b) Thermal output of CHPs in comparison to thermal load. 
c) Electric output of CHPs, electric load and interactions with electricity market. 
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Figure 6.13: Hourly averaged results for 2 CHPs scenario with hybrid discrete-continuous time model. 
a) Supply temperature at CHP1. b) Thermal output of CHPs in comparison to thermal load. 
c) Electric output of CHPs, electric load and interactions with electricity market. 
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In comparison to a simulation using the node method (Section 2.1.2) with the opti- 
mization results of supply temperature of CHP1, thermal output of CHP2 and thermal 
load as input, the thermal output of CHP1 in the optimization result has a root mean 
square deviation of 32.850. 


6.1.4 Optimization results with delay matrix approach 


In the following, the case study results of an optimization of the heating grids de- 
scribed in Section 6.1.1 with the delay matrix approach introduced in Chapter 5 are 
presented. As in the previous section, thermal losses are not considered being of 
minor importance for the studied cases. 


The computations for the delay matrix approach were run on the setup used in the 
previous section (an Intel Silver 4110 CPU with 2.10 GHz and 16 GB RAM using 
Julia/JuMP [Ran16] for modeling and Gurobi [Gur16] as MILP solver). 


One CHP case 


Figure 6.14 shows the results of the delay matrix approach for the one CHP case being 
achieved in 0.035 s leading to an overall cost of 291,173.4 €. 


The electric output of the CHP presented in Figure 6.14 c) alternates between high 
and low electric output in the early morning hours until hour 6. As for the previous 
optimization results this is caused by an electricity market price swinging around the 
marginal production cost. In Figure 6.14 b) there is a thermal discharge (CHP output 
below load) of the grid around hours 4, 8, 12 and 18, being the high price hours. 
The grid is charged (CHP output above load) around hours 2, 6, 11 and 13, where 
electricity prices are lower. This charging and discharging is reasonable as a lower 
thermal output allows the CHP to increase electric output (see Figure 2.1) and, thus, 
revenues are increased as electric production is moved to higher price times. This load 
shifting can as well be observed in Figure 6.14 a) showing the supply temperature 
at the CHP. The supply temperature is increased to charge the grid and decreased 
to discharge the grid. At the end of the time horizon after the last price peaks, the 
supply temperature remains at the lowest level, discharging the thermal energy stored 
in the heating grid completely. This is reasonable as the model formulation does not 
account for thermal energy stored in the heating grid at the end of the time horizon. 


A root mean square deviation of 49.159 for the thermal output of the CHP is obtained 
when comparing the optimization result to a simulation that uses the node method 
(Section 2.1.2). The inputs of this simulation are the optimization result of the supply 
temperature of the CHP and the thermal load profile. 
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Figure 6.14: Results for 1 CHP scenario with delay matrix approach. a) Supply temperature at CHP. 
b) Thermal output of CHP in comparison to thermal load. c) Electric output of CHP, elec- 
tric load and interactions with electricity market. 
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Two CHPs case 


The optimization result for the two CHPs scenario as introduced in Section 6.1.1 using 
the delay matrix approach (Section 5) is shown in Figure 6.15. This result with an 
objective value of 349,537 € was achieved in 111 s having a MILP gap of 0.0044 %. 


Comparing Figure 6.15 with the results for the one CHP scenario in Figure 6.14 we 
observe a similar pattern of the supply temperature: Before hours 8, 12 and 18 the 
supply temperature of CHP1 is increased to store thermal energy in the heating grid, 
whereas at times 8 to 10, 12 and starting in hour 18 the supply temperature of CHP1 
drops back to its original level, discharging the heating grid. In these periods of 
discharging thermal energy from the heating grid the thermal output of CHP1 can 
be reduced and its electric output can thus be increased (see Figure 2.1). The overall 
benefit is increased by this charging pattern increasing the electricity output at higher 
price times. For CHP2 there is a linear connection between electric and heat output. 
Thus, it runs at all times when it is profitable to sell electricity to the market, except 
in hour 11 where as much energy as possible should be stored in the heating grid by 
CHP1. In hour 7 the electricity price jumps above the cost of production of CHP1 and, 
thus, it increases electric output to its maximum level. 


Running a simulation of the heating grid case using the node method (Section 2.1.2) 
with optimization results of supply temperature of CHP1, thermal output of CHP2 
and thermal load as input, results in a root mean square deviation of 51.160 for the 
thermal generation of CHP1 between simulation and optimization result. 


6.1.5 Reference case: sequential approach 


As a reference case we use the approach proposed in [SLM*17] to compare its per- 
formance with the approaches of this thesis. In this work, the optimization problem 
is split into, first, a MILP unit commitment problem that does not consider thermal 
dynamics of the heating grid and, second, a NLP economic dispatch problem with 
fixed unit commitment but considering the thermal dynamics of the grid [SLM*17]. 
For the unit commitment problem the formulation presented in Section 2.1.1 is used. 
As our cost function is linear without quadratic terms we get a MILP formulation 
of the unit commitment problem in comparison to [SLM*17] using a quadratic cost 
function. For the economic dispatch problem with fixed unit commitment the primal 
problems of the global optimization presented in Section 3.2 are used. Like for the 
global optimization, thermal losses are considered in the NLPs. 


The reference case is run on an Intel Xeon CPU E5-2650 v3 with 2.3 GHz, 16 cores and 
16 GB RAM using the modeling environment Julia/JuMP [Ran16] with Gurobi 8.1 
[Gur16] as solver for the MILPs and IPOPT [WBO06] as local solver for the NLPs. 
This setup was used for the global optimization case study, too. The results for the 
sequential approach are achieved in at most 30 min. Major solution time is required 
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Figure 6.15: Results for 2 CHPs scenario with delay matrix approach. a) Supply temperature at CHP1. 
b) Thermal output of CHPs in comparison to thermal load. c) Electric output of CHPs, electric 
load and interactions with electricity market. 
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for the NLPs. This is significantly longer than the NLP solver runs in the global 
optimization scheme due to worse initial values. In the global optimization scheme, 
the initial values are taken from the results of the dual solution and, hence, are closer 
to the solutions than with standard initial values. This case study for the reference 
case is published in [MLH19]. 


One CHP case 


If only one CHP is available, it needs to run throughout the full optimization horizon 
as the CHP is the only source to supply the heat demand. Thus, the CHP is always 
running and solving a unit commitment optimization is obsolete. The economic dis- 
patch problem is run directly without any special initial conditions for IPOPT [WB06]. 
The results for the optimization runs for the finite-volumes model (Model 3.2) as well 
as the hybrid pipe model (Section 3.3) approach the same solutions as the correspond- 
ing global optimization scheme using those models as primal problems. 


Accordingly, Figure 6.6 and Figure 6.7 showing the result of the global optimization, 
display the results of the reference case using the corresponding primal problem as 
economic dispatch problem, too. Running only the economic dispatch NLP takes 
about 2 to 3 seconds. 


Two CHPs case 


For the scenario with two CHPs the sequential approach following [SLM*17] does 
not yield a feasible solution. The MILP unit commitment problem solves without 
any issues, but when running the economic dispatch NLPs, the local solver (IPOPT 
[WB06]) always converges to an infeasible solution for both NLP formulations, the 
finite-volumes model (Model 3.2) as well as the hybrid pipe model (Model 3.3). This 
is most likely caused by the start values of IPOPT: As the unit commitment problem 
does not compute any suitable starting values for the economic dispatch problem be- 
sides the commitment status, standard start values are taken. In contrast, the primal 
problems in the global optimization are initialized with the results of the dual prob- 
lem. Hence, this issue does not arise in the cases presented in Section 6.1.2 which use 
multiparametric delay modeling with the same non-linear model formulations. 


This situation might improve with an algorithm finding better start values. However, 
the result will differ from the results of the global optimization, as the commitment 
status of CHP2 after the unit commitment optimization in the sequential approach is 
different from the commitment status of CHP2 achieved in the global optimization 
runs as shown in Figure 6.16. 


To further compare the sequential approach with the global optimization algorithm 
an additional case study with both CHPs at the location of CHP1 is run. The other pa- 
rameters of the two CHP scenario remain. Thus, CHP1 and CHP2 are jointly feeding 
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Figure 6.16: Comparison of unit commitment of CHP2 between reference case and global optimization 
approaches 


into the 10 km pipe, which transports the thermal energy to the load area. Here, for 
the global optimization with both primal models, we get the same unit commitment 
decision as with the unit commitment in the sequential approach. In addition, with 
the finite-volumes model in the sequential approach the economic dispatch converges 
to the same overall solution as the global optimization scheme with the finite-volumes 
model as primal problem. However, when using the hybrid pipe model for economic 
dispatch in the sequential approach and as primal problem in the global optimization 
scheme, the algorithms result in different solutions. As the NLPs are only solved with 
a local solver, this is caused by the different starting values for the optimization of the 
hybrid pipe model in the sequential scheme and in the global optimization scheme. 


6.1.6 Discussion of results and comparison of modeling approaches 


Overall, the optimization results of the different approaches show several similarities: 
Electric generation of the CHP turns profitable in hour 6 in the one CHP case such that 
its electric output is increased to the maximum possible value for the current thermal 
output. In the two CHPs case the electricity generation becomes profitable one hour 
later in hour 7. In the one CHP scenarios there is a large fluctuation of the electric 
generation in the early morning hours as the electricity price is swinging around the 
marginal cost independent of the optimization approach. 


Furthermore, in all results, the thermal output of the coal fired power plant is reduced 
around high price hours (hour 8 and hour 18) to allow a higher output of electric en- 
ergy. Here, the thermal load is supplied by thermal energy stored in the heating grid 
by an increase of supply temperature in previous time steps. The supply temperature 
is decreased to its minimum level in order to discharge the thermal energy stored in 
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the heating grid. The gas fired power plant in contrast is run in all higher price hours 


at maximum capacity as electric and heat output are directly linked. 


With all methods, the optimizer pushes the supply temperature at the CHP to its 
boundaries leading to large supply temperature changes. In real-world situations 
large variations in temperature are usually not wanted as larger variations in temper- 
ature cause higher changes in volume requiring bigger buffer tanks and can increase 


material fatigue. 


Table 6.2: Results overview for one CHP case 


Delay matrix Hybrid time Global: Global: 
grid finite- hybrid pipe 
volumes 
Temperatures in °C 60 to 130 90, 95, 100, 60 to 130 60 to 130 
(Continu- 105, 110, (Continuous) (Continu- 
ous) 115, 120, 130 ous) 
Objective in € 291,173.4 292,565.9 291,235.7 294,386.0 
Lower Bound in € - - 288,371.3 288,265.1 
(global) 
Solution Time in s 0.035 800 146.19 125.00 
MILP Gap in % 0 0.0513 - - 
Global Gap in % = z 1.0 2.1 
Number of Vari- 403 Cont. 3,481 Cont. 8,738 Cont. 19,764 Cont. 
ables (MILP) 124 Int. 1,319 Int. 11,026 Int. 11,026 Int. 
Table 6.3: Results overview for two CHPs case 
Delay matrix Hybrid time Global: Global: 
grid finite- hybrid pipe 
volumes 
Temperatures in°C 60 to 130 90, 95, 100, 60 to 130 60 to 130 
(Continu- 105, 110, (Continuous) (Continu- 
ous) 115, 120, 130 ous) 
Objective in € 349,537 346,490 348,252.9 351,112.0 
Lower Bound in € - - 330,871.1 334,785.9 
(global) 
Solution Time in s 111 800 1,637.0 2,907.1 
MILP Gap in % 0.0044 0.0554 - - 
Global Gap in % - - 4.99 4.6 
Number of Vari- 506 Cont. 3,581 Cont. 9,489 Cont. 11,337 Cont. 
ables (MILP) 447 Int. 1,462 Int. 11,676 Int. 15,076 Int. 
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Comparing the results of the different methods for both cases studied in Table 6.2 and 
Table 6.3, we can observe that the objective values have the same order of magnitude 
and vary comparably little. This shows, that thermal losses, considered in the global 
optimization approaches, but not considered in the delay matrix and the hybrid time 
grid approach, are of minor importance for the studied scenarios with one pipe with 
a large diameter. The lower bounds achieved with the global optimization approaches 
are comparable, too. For the one CHP case, the global optimization approaches show 
very good performance. For the two CHPs case, solving the dual problem becomes 
a lot more challenging such that solution time and optimality gap increase. Both 
however stay at a reasonable level. 


Assuming delay times are given, the delay matrix approach results in the smallest 
optimization model offering the fastest solution times. However, the extension to a 
second CHP presented in Section 5.3.2 increases the solution time by multiple orders 
of magnitude as it introduces many binary variables. Nevertheless, it still remains the 
fastest solver for the two CHPs case. Problem size as well as solution times for the 
hybrid time grid model are larger but still within an acceptable range for day ahead 
planning problems. As this approach offers an exact modeling of the transport delays, 
it produces results closer to the real-world situation as the delay matrix approach. 
However, possible temperatures are limited to a discrete set of temperature levels. 
For the case studies, steps of 5 K are assumed which should be accurate enough 
for most real-world applications. This was validated with a simulation model: To 
evaluate solution accuracy the optimization result of the thermal output of CHP1 was 
compared with a detailed simulation using the optimized supply temperature profile 
as input. The root mean square deviation between optimization result and detailed 
simulation reduces from 49.159 to 29.431 for the one CHP case and from 51.160 to 
32.850 for the two CHPs case, if the hybrid time grid model is used instead of the 
delay matrix approach. 


The temperature dynamics in the primal solutions using the finite-volumes model in 
the global optimization scheme are not completely realistic, as the outflow tempera- 
ture of the pipe varies less than the inflow temperature due to the smoothing effects 
of the chosen models. The hybrid pipe model leads to a more accurate representa- 
tion of the transport delay and to a larger objective value as shown in Table 6.2 and 
Table 6.3. This is reasonable as a higher variability of the CHP supply temperature 
allows to reach a better objective value of the primal problems with the finite-volumes 
model. 


The sizes of the dual problems in the global optimization algorithm strongly depend 
on the number of iterations the algorithm needs to achieve a good solution. Solution 
times for the one CHP case are competitive in comparison to the hybrid time grid 
approach. The solution gap varies from 1.0 % to 2.1 % depending on the chosen 
primal problem being achieved in less than 2.5 min. However, for the two CHPs case 
solution times increase significantly even though the accepted global optimality gap 
is increased to 5 %. About 30 min and 50 min are needed to reach this gap with the 
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different primal problems. If the MILP solver time limit for solving the lower bound 
problem is increased, the global optimality gap can be closed further. Nevertheless, 
starting from iteration 3 the MILP solver struggles to find a feasible solution for the 
lower bound problem. Thus, warm-starting the solving of the lower bound problem 
might improve this situation. Especially when larger heating grids shall be optimized, 
such improvements to the algorithm are required. 


A unit commitment optimization using the basic MILP model presented in Section 
2.1.1 not considering neither thermal dynamics of the heating grid nor thermal losses 
reaches an objective value of 297,788.6€ for the one CHP case and an objective value 
of 355,059.0€ for the two CHPs case. Thus, heating grid operators could achieve cost 
savings between about 3,500 € and 6,700 € per day translating to relative savings of 
about 1.2 % to 3.3 % using the heating grid as thermal storage. 


6.2 Real-world case study 


One of the goals of this thesis is to find optimization approaches for district heating 
grids which can directly be applied to real-world cases. As the delay matrix approach 
introduced in Chapter 5 enables a fast and efficient optimization of large grids, it is 
applied to the district heating grid of the city of Kiel, Germany in this section. This 
real-world case study is published in [MKV*19]. 


All calculations of this real-world case study are run on an Intel Xeon v2 with 2.6 GHz 
(2 CPUs) and 4 GB RAM using the modeling environment Julia/JuMP [Ran16] and 
MILP solver Gurobi [Gur16]. 


6.2.1 The district heating grid of Kiel 


Kiel is a city at the German coast of the Baltic Sea with about 250,000 inhabitants. 
The district heating grid of Kiel has a total pipe length of about 370 km serving about 
7,000 building connection points [Stal5b]. With this size it is a very representative 
large heating grid (defined as grids above 100 km length) being only slightly longer 
than the average large heating grid in Germany [SDEW12]. The largest heat supplier 
in Kiel was a coal fired power plant with about 354 MW electric output which was 
replaced by a gas engine plant with twenty gas combustion engines each having an 
electric and heat output of about 10 MW in 2019 [Stal5b, Sta16]. 


For this case study, 37 days in the periods January to March and November to Decem- 
ber are selected from a full year of measurement data as in [MKV‘19]. These days 
have been chosen, as here the coal fired power plant supplies almost all heat demand 
in the heating grid and, thus, flow directions are comparably stable. Therefore, it is 
possible to identify a stable delay matrix. Based on measurement availability, several 
consumption zones are introduced which split the demand of the heating grid among 
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them. For each consumption zone, its share of the overall consumption as well as the 
transport delay from the generation site are calculated allowing to create the delay 
matrix MI], t] introduced in Table 5.1 of Chapter 5. The transport time from genera- 
tion site to the main consumption area including the city center is about one to two 
hours for the studied days. Usually, all consumption areas in the heating grid are 
reached by the hot water after at most 8 hours. 


The formulation introduced in Section 5.3.1 is used to consider the additional heat 
losses caused by the supply temperature increases chosen in the optimization. The 
heat loss factor &joss is calculated using topology information giving installed lengths 
and diameters of all pipes in the heating grid in combination with data sheet in- 
formation on the heat loss coefficient of the pipes. For some (mostly older) pipes, 
which have been insulated during the installation works, no data sheet information 
is available. Thus, the heat loss coefficient largely varies and is worse than the heat 
loss coefficient of pre-insulated pipes. Based on consultations with Stadtwerke Kiel, 
the heat loss coefficient of those pipes is assumed to be the double of the heat loss 
coefficient of a pre-insulated pipe with the same diameter. 


For every day considered, prices from the German/ Austrian price zone of the EPEX 
spot day ahead market [EPE] are taken as electricity prices for sales and with a minor 
offset for purchases. In the optimization these hourly varying prices are considered 
as parameters. 


Two different generation scenarios are used, which are inspired by past and new main 
generators in the heating grid of Kiel, but do not reflect the exact real-world size and 
setting. As in both scenarios only one CHP supplies the full demand, the thermal 
output of the original CHPs is increased by adjusting size or number of units. In the 
first scenario, the district heating grid is supplied by a coal fired power plant with an 
extraction condensing turbine, 130 °C maximum supply temperature and maximum 
output of 300 MW electric and 300 MW thermal power. In part-load operation, only a 
minor drop of efficiency is assumed. As an extraction condensing turbine is installed, 
electric and thermal output are linked with the feasible operation region shown in 
Figure 2.1. 


In the second generation scenario the district heating grid is supplied by a gas engine 
plant consisting of thirty gas combustion engines each having a maximum thermal 
output of 10 MW and a maximum electric output of 10 MW. The gas combustion 
engines supply heat at a maximum temperature of 115 °C having a linear power- 
to-heat ratio. In part-load operation, a linear efficiency drop is assumed, having the 
highest efficiency at maximum output. The start up cost of the gas engines is assumed 
to be low. An overview on the different parameters of the generation scenarios used 
in the real-world case study is given in Table 6.4. 
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Table 6.4: Generation scenarios in real-world case study 


Coal fired plant Gas engine plant 


Number of units 1 30 
Electric output per unit in MW < 300 < 10 
Thermal output per unit in MW < 300 < 10 
Supply temperature in °C < 130 <115 
Part-load efficiency drop small big 
Power to heat ratio see Figure 2.1 constant 


6.2.2 Results for a December day 


To show detailed results for one day, the delay matrix approach introduced in Section 
5 is run for a December day for both generation scenarios presented in the previous 
section. Using measurement data of the district heating grid of Kiel from this date, 
the overall thermal energy consumption profile as well as the share of consumption 
of the different consumption zones and the time delay from thermal generation to 
these consumption zones are calculated. This information is used to parameterize the 
delay matrix M|}, t]. All electric energy produced by the CHPs is sold at the energy 
market. Figure 6.17 shows the relative price profile for this December day scaled with 
the price peak of the day. 


1 


0.8 


0.6 


0.4 


0.2 - | 


Relative electricity price 


| | 
2 4 6 8 10 12 14 16 18 20 22 24 
time of day inh 


Figure 6.17: Relative electricity price for a December day 


It is assumed that the supply temperature before optimization is at the minimum 
possible temperature according to the technical connection conditions of the heating 
grid. Thus, the supply temperature profile before optimization is calculated based on 
the outdoor air temperature. For the second generation scenario with the gas engine 
plant it is additionally assumed that a temperature increase cannot exceed 10 K. 
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Figure 6.18: Result of real-word case study for coal fired power plant (one day). a) Supply temperature 
increases. b) Thermal energy stored in the district heating grid with supply temperature vari- 
ations and additional thermal losses caused by supply temperature increase. 


The optimization result for the first generation scenario consisting of the coal fired 
power plant with extraction condensing turbine is shown in Figure 6.18. Thermal 
energy is stored in the heating grid at hour 2, hours 6 to 7 and hours 15 to 16 with 
an increase of supply temperature at the producer. These periods are followed by a 
discharge of thermal energy from the heating grid. During this discharging phase, the 
thermal output of the coal fired CHP is reduced enabling a higher output of electricity 
of the CHP according to the feasible operation region of this CHP (Figure 2.1). Thus, 
the supply temperature increases to charge the heating grid before peaks of electricity 
price around hours 8 and 18. During these peak times the thermal energy stored 
within the heating grid is used to supply the thermal demand and the CHP can 
increase its electric output. Hence, electric generation of the CHP is moved to hours 
where higher revenues are possible. The increase of thermal losses is comparably low. 
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Nevertheless, it reduces overall profits of storing thermal energy in the heating grid. 
This optimization using the delay matrix model finished in less than 5 ms increasing 
the overall profit of operations by several hundred Euros (~0.2 %) in comparison to a 
scenario without thermal storage. 
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Figure 6.19: Result of real-word case study for gas engine plant (one day). a) Supply temperature increases. 
b) Thermal energy stored in the district heating grid with supply temperature variations and 
additional thermal losses caused by supply temperature increase. 


The result for the second generation scenario consisting of a gas engine plant with 30 
individual gas combustion engines is shown in Figure 6.19. In contrast to the coal fired 
plant, output of thermal energy increases with output of electricity for the gas engine 
plant. Hence, the optimizer increases supply temperature at high price times around 
hours 8 and 18 storing thermal energy in the grid and increasing power output of the 
CHP. However, there are many other variations in supply temperature which cannot 
be explained with the volatile energy price. They are caused by the drop of efficiency 
in part-load operations. As it is more efficient to run a gas combustion engine at 
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maximum speed, the optimizer pushes to run gas units at maximum output if they 
are running. But as the thermal demand is usually not at an integer multiplier of 
the gas engine size, the thermal storage capacity of the heating grid is used to buffer 
the differences. Thermal losses are less than for the coal fired power plant case as 
the temperature increases are smaller. Despite these rather small supply temperature 
increases, profits for the gas engine plant case are higher than for the coal fired power 
plant case, achieving savings of several thousand Euros in this day (~2.5 % of overall 
generation cost). Looking at the end of the optimization horizon it must be noted 
that there is a rather late big supply temperature increase in hour 22, such that the 
heating grid is not fully discharged at the end of the optimization horizon. Thus, if 
this thermal energy could be used e.g. extending the optimization horizon, savings 
might increase further. As more integer variables are involved, the optimization of 
the gas engine plant case takes 1.7 s which is several orders of magnitude longer than 
the optimization of the coal fired power plant case. Nevertheless, this remains an 
acceptable solution time. 


6.2.3 Evaluation of economic potential 


To evaluate the economic potential of using the heating grid as thermal storage, the 37 
days with feasible measurement data of the heating grid of the city of Kiel are grouped 
into several scenarios. Grouping by different optimization horizons we get 


e 37 scenarios for a one day (24 h) horizon, 
e 27 scenarios for a two day (48 h) horizon, 
e 19 scenarios for a three day (72 h) horizon and 
e 14 scenarios for a four day (96 h) horizon. 


Note that the scenarios are partly overlapping for multi day scenarios. Several cases 
are studied running optimizations for both generation setups presented in Section 6.2.1 
and for all scenarios mentioned above. 


For the coal fired power plant the following four cases are run with differences in 
the supply temperature level assumed before optimization and in consideration of 
thermal losses. 


e minimum supply temperature without thermal losses 
e minimum supply temperature with thermal losses 
e measured supply temperature without thermal losses 


e measured supply temperature with thermal losses 
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The minimum supply temperature is calculated with the outdoor air temperature ac- 
cording to the technical connection conditions of the heating grid of Kiel [Stal5a]. 
Measured supply temperature cases take the actual measurement of supply tempera- 
ture at the CHP in the scenario as optimization input. 


For the gas engine plant three cases are studied. All of them consider thermal losses 
and are using the minimum supply temperature calculated with the outdoor air tem- 
perature according to the technical connection conditions of the heating grid of Kiel 
[Stal5a]. It is not feasible to use the measured temperature profiles as they often are 
above the maximum supply temperature of the gas engine plant. The cases studied 
for the gas engine plant are: 


e original part-load efficiency 
e higher part-load efficiency 


e original part-load efficiency and limitation of supply temperature increase to 10 K 


Coal fired power plant 


An overview on the case study results for the coal fired power plant is given in Fig- 
ure 6.20. Comparing these results, three major influences on the relative savings of 
using the heating grid as thermal storage in comparison to a scenario without thermal 
storage can be identified. 


First, it can be noted that increasing the optimization horizon increases the average 
and decreases the variance of the savings for all cases. With a longer time horizon 
the variance in the electricity prices between the scenarios decrease and, accordingly, 
the variance of the relative savings decreases: If there is one very favorable electricity 
price increase or decrease in a day, its impact is reduced, if not only this day, but 
also the neighboring days without such a favorable electricity price development are 
considered. The increase of the average and median relative savings with increasing 
time horizon is most likely caused by inefficient usage of the thermal storage capa- 
bility of the heating grid at the end of the optimization horizon: If the heating grid 
is not fully discharged at the end of the optimization horizon, the remaining thermal 
energy is considered lost in the optimization. Thus, storing energy at the end of the 
optimization horizon is less interesting and changes in electricity price at the end of 
the optimization horizon cannot be used efficiently. Having an optimization horizon 
spanning multiple days allows to leverage the variations of electricity price in the 
evening for all but the last day without limitations. Thus, the negative impact of this 
problem decreases. 


Second, it can be observed that the cases with minimum supply temperature are 
creating higher benefits than the cases with measured supply temperature. As the 
measured supply temperature always lies above the minimum supply temperature, 
the maximum increase of supply temperature before reaching the technical limitations 
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Figure 6.20: Relative savings achieved in real-world case study for the coal fired power plant for different 
time horizons and cases: minimum supply temperature without thermal losses (blue) and with 
thermal losses (gray); measured supply temperature without thermal losses (orange) and with 
thermal losses (yellow). The boxes mark the lower and upper quartiles, hence 50 % of the data 
lies inside the box. The crosses mark the average and the line inside the box the median of 
the data. The whiskers show the variability of data within at most a distance from lower and 
upper quartile of 1.5 times the inner quartile range. Points denote outliers being data outside 
this range. 


of the coal fired power plant is smaller. The thermal energy storage capacity of the 
heating grid is proportional to the possible increase of supply temperature. Thus, 
having a higher supply temperature before the optimization, reduces the thermal 
storage capacity of the heating grid and, thus, less savings can be achieved. 


Third, considering heat losses reduces the economic benefits. This connection is ex- 
pected, as considering heat losses requires additional thermal generation increasing 
cost without increasing revenues. As can be seen in the one day example in Sec- 
tion 6.2.2, temperature increases in the coal fired power plant case can be high (up 
to 30 K). Thus, the increase of thermal losses caused by this temperature increase re- 
duces the relative savings. Thermal losses induced by temperature increases play an 
important role for the considered heating grid in the city of Kiel, as the temperature 
increases reach all pipes in the system, even the very small ones connecting single 
houses. Hence, the surface of the heating grid network is very large in comparison to 
its volume. For heating grids with only larger pipes the importance of thermal losses 
is lower. The same applies to heating grids where the temperature in large pipes and 
small pipes is decoupled. Here the storing of thermal energy can be limited to the 
larger pipes where temperature increases lead to smaller thermal losses than in small 
pipes with high surface to volume ratio. 
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The solution times for all scenarios with the coal fired power plant are extremely fast, 
never exceeding 50 ms. 


Gas engine plant 


Figure 6.21 shows an overview on achievable relative savings based on different sce- 
narios and cases for the gas engine plant. Only cases considering thermal losses and 
assuming the minimum supply temperature as pre-optimization supply temperature 
are calculated for the gas engine plant. A similar trend as for the coal fired generation 
can be observed regarding an increase of time horizon: In scenarios with a longer time 
horizon, the average and median achieved relative savings increase and the variance 
of the relative savings decrease. However, the second effect is less important, as for 
the gas engine plant major savings are achieved running the combustion engines at 
maximum efficiency rather than shifting generation based an varying energy prices 
as discussed for the one day scenario in Section 6.2.2. 
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Figure 6.21: Relative savings achieved in real-world case study for gas engine plant for different time hori- 
zons and cases: original part-load efficiency (light blue); higher part-load efficiency (green); 
original part-load efficiency and maximum temperature increase of 10 K (red). The boxes 
mark the lower and upper quartiles, hence 50 % of the data lies inside the box. The crosses 
mark the average and the line inside the box the median of the data. The whiskers show the 
variability of data within at most a distance from lower and upper quartile of 1.5 times the 
inner quartile range. Points denote outliers being data outside this range. 


This implication, that the main benefits for the gas engine plant scenarios are drawn 
from running the individual gas units at maximum efficiency, is supported if the cases 
with original part-load efficiency and higher part-load efficiency are compared. The 
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achieved benefits for the gas engine plant with higher part-load efficiency are signif- 
icantly smaller than the relative savings achieved with original part-load efficiency. 
Thus, for the gas engine plant cases using all running units at maximum efficiency is 
the major reason for the savings generated. 


In the third case studied for the gas engine plant, the temperature increase was lim- 
ited to 10 K. Here, the achieved relative savings are only slightly reduced. Thus, the 
temperature increases required to achieve benefits with the gas engine plant are com- 
parably small and a lot less important than the temperature increases calculated for 
the coal fired power plant. This is supported by the results for one example day dis- 
cussed in Section 6.2.2. Despite this fact, the savings generated using the heating grid 
as a thermal storage with the gas engine plant are importantly higher than the savings 
generated with the coal fired power plant shown in Figure 6.20. This is remarkable as 
thermal losses are considered for all gas engine scenarios. The major reason for this 
is most likely that the coal fired plant can generate benefits only by shifting genera- 
tion between different price time slots, which requires a rather large thermal storage. 
The gas engine CHP can additionally generate savings by buffering thermal energy to 
run its units at maximum efficiency if they run, which requires less thermal storage. 
The gas engine plant has also a more direct link between thermal and electric output 
than the coal fired CHP. Hence, changes in thermal output of the assumed gas engine 
plant have stronger influence on its electric output than for the coal fired power plant. 
Thus, smaller changes in thermal generation are required to achieve the same change 
in electric output. 


Note that part-load efficiencies and generation cost of the gas engine plant and the 
coal fired power plant are estimated. Real-world behavior of a gas engine plant or a 
coal fired power plant might vary. 


Finally, in Figure 6.21 it can be noted that for the 24 h and 48 h scenarios there 
are always some scenarios where no benefits can be achieved at all. This can be 
explained by Figure 6.22 which shows that the days not achieving any benefits are 
days with very low ambient air temperatures. At low ambient air temperatures the 
heating grid must be supplied with high supply temperatures (at least 115 °C for 
ambient air temperatures below -5 °C according to the technical connection conditions 
[Stal5a]). Hence, there is no room for a temperature increase as the gas engine plant 
can only produce up to 115 °C and no thermal energy can be stored in the heating 
grid. Thus, the energy storage capabilities of the heating grid can mainly achieve 
interesting economic savings at average ambient air temperatures above 0 °C. 


Solution times for the scenarios with the gas engine plant are slower than optimiza- 
tions for the scenarios with the coal fired power plant, as more integer variables rep- 
resenting the commitment status of the individual gas engines are involved. For the 
scenarios and cases presented in this section they are between 1 s for the 24 h scenar- 
ios and reach up to 100 s being the MILP solver time limit for some of the 72 h and 
96 h scenarios. 
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Figure 6.22: Dependency of relative savings for the gas engine plant on daily average 
ambient air temperatures for scenarios with 24 h horizon 


6.3 Summary of chapter 


This chapter contains the evaluation of the methods developed in this thesis. In Sec- 
tion 6.1 all optimization methods developed in this thesis were compared using two 
small test scenarios: one with one CHP and one with two CHPs. The delay matrix 
approach (Chapter 5) proves its competitiveness being the fastest optimization model 
in all studied cases, very often providing a solution in less than a second. Solving 
the hybrid time grid model (Chapter 4) takes longer, but solution times stay within 
an acceptable range for a day ahead planning problem being at most 10 minutes. 
Comparing the optimization result with a simulation using the optimization result as 
input, the results of the hybrid time grid model show a higher precision than the re- 
sults using the delay matrix approach. This is expected as the hybrid time grid model 
allows an accurate representation of the time delay. 


The global optimization approach (Chapter 3) was run using both introduced primal 
problems. Solution times of the global optimization approach are competitive only 
for the one CHP case. For the two CHPs case, solution times increase to 30 min and 
50 min. The same applies to the global optimality gap achieved: for the one CHP 
case gaps of 1 % and 2.1 % are possible, whereas for the two CHPs case gaps of about 
5 % are achieved. Thus, increasing the size of the studied heating grid increases the 
computational effort drastically. The objective value with the finite-volumes model as 
primal problem is always slightly smaller than the objective value with the hybrid pipe 
model as primal problem. This can be explained by the less accurate representation 
of the time delay in the finite-volumes model: As there is a smoothing of temperature 
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introduced by the finite-volumes, the temperature at the CHP can vary more than with 
the hybrid pipe model (sometimes even to values not possible in real-world). Thus, 
variations in generation are larger allowing to create more benefits. In addition to the 
approaches presented in this thesis, the sequential approach proposed in [SLM*17] 
was run for both studied cases. This method reaches the same result as the global 
optimization approaches for the one CHP case, but does not succeed in delivering 
a solution for the two CHPs case. Overall, all approaches give comparable results, 
being slightly above the calculated global lower bounds. 


In Section 6.2 the delay matrix approach (Chapter 5) was applied in a real-world case 
study using measurement data from the heating grid of the city of Kiel. The delay 
matrix was created using transport delays to and share of consumption of load areas, 
being extracted from this measurement data. A thermal loss factor was identified 
using grid topology and pipe data sheet information. Two different generation cases 
were considered: a coal fired power plant with extraction condensing turbine and a 
gas engine plant consisting of 30 gas engines. Different scenarios were used based 
on the real-world measurement data of 37 days in the periods January to March and 
November to December. All results were compared with an optimization not consid- 
ering the storage capabilities of the heating grid and corresponding relative savings 
were calculated. After the presentation of results for one day in December in Sec- 
tion 6.2.2, a summary of all studied scenarios and cases was presented in Section 6.2.3. 
Several influences on the achievable savings using the heating grid as a storage have 
been identified: 


e A longer optimization horizon allows higher savings, as inefficiencies at the end of 
the optimization horizon occur less often. 


e Ambient air temperature has an important influence on achievable savings, because 
no increases in supply temperature are possible when the heating grid operates 
close to full load. 


e Heat losses decrease the achievable savings. If smaller pipes are decoupled from 
large pipes (with heat exchangers), this effect can be reduced. 


e Shifting load/generation between time slots requires a large amount of thermal 
storage, hence, temperature variations and the influence of thermal losses are big. 


e Using the storage capacity to run plants at maximum efficiency requires less ther- 
mal storage. Thus, temperature variations are smaller, thermal losses are less im- 
portant and savings are higher. 


All in all, the delay matrix approach proved to be capable to optimize real-world 
heating grids. Optimization times stayed in an acceptable range mostly being solved 
in less than a second. Only few larger scenarios with longer optimization horizon 
took up to 100 s to solve. The resulting savings are up to 6 % or a few thousand 
Euros per day and, thus, could easily reach more than 180,000 € per winter season 
(assuming average savings of 1,000€ per day for half a year). However, they are 


6.3 Summary of chapter 103 


dependent on the grid topology as well as on the generation setup. Hence, no general 
recommendation on using the heating grid as thermal energy storage can be given 
and an analysis of generation setup and grid topology is required to identify saving 
potentials. 


7 Conclusions 


With increasing share of renewable power generation in electric grids, conventional 
generation needs to become more and more flexible to react to fluctuating renewable 
feed in and volatile electricity prices. As heat demand in heating grids is usually not 
flexible, CHPs are facing major challenges. One way to enable a more flexible oper- 
ation of CHPs is the installation of thermal storage tanks. A promising alternative, 
that avoids installation effort and investment cost, could be to directly use the abil- 
ity of the thermal dynamics of the heating grid to store thermal energy. However, 
considering the thermal dynamics of a heating grid in operations planning of CHPs 
is complex due to bilinear terms and a variable-dependent time delay. Several opti- 
mization algorithms in literature address this non-convex optimization problem, but 
no approach reaches the global optimum. In addition, such optimization of energy 
flows considering heating grid dynamics is not widely used in practice. 


Multiparametric delay modeling presented in Chapter 3 is the first optimization al- 
gorithm proving global optimality of its solution. It introduces a new outer approxi- 
mation of the variable-dependent time delay which is formulated as MILP problem. 
Thus, it can be solved with off-the-shelf solvers. In combination with models of ther- 
mal dynamics from literature sources, this outer approximation is used in a primal- 
dual optimization scheme. Using this iterative approach, the solution converges to the 
global optimum. In the case studies in Chapter 6 it is shown that for a small heating 
grid with one CHP and one load the algorithm finds a solution with a global optimal- 
ity gap of 1 % in less than two minutes. However, as soon as problem size increases 
(e.g. adding one CHP) solution times increase importantly. Thus, multiparametric de- 
lay modeling is a great approach to benchmark other optimization approaches with 
respect to global optimality, but without further improvements it is not a suitable 
solution for every day operations planning of a heating grid. 


To reach a single, very exact model of the thermal dynamics in a heating grid, which 
can be used for every day operations planning, the hybrid time grid approach was pre- 
sented in Chapter 4. It adapts optimization models originally developed for pipeline 
scheduling in the petrochemical industry. To suite the needs of heating grid optimiza- 
tion a hybrid discrete-continuous time grid replaces the continuous time grid of the 
pipeline scheduling model and discrete temperature levels for the supply tempera- 
ture are introduced. This of course limits the choice of supply temperatures, but if 
the number of temperature levels is high enough (e.g. every 5 Kelvin) the influence 
of this discretization is limited. With these discrete temperature levels and the hybrid 
discrete-continuous time grid, a MILP model formulation of the thermal dynamics 
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is possible. The resulting MILP optimization model for the full problem setup can 
be solved directly with standard solvers. Solutions times in the case studies (Sec- 
tion 6.1.3) are acceptable. An extension to heating grids with several interconnected 
pipes and reversible flow is possible, as extensions of the underlying scheduling mod- 
els for pipeline scheduling already support such pipe networks. However, all pipe pa- 
rameters need to be known to correctly parameterize the grid. Thus, the hybrid time 
grid approach seems to be very suitable for transport grids with a limited number of 
pipes. 

For distribution grids with a high number of pipes with often unknown parameters, 
a different approach is needed. In comparison to most literature approaches, the 
delay matrix approach presented in Chapter 5 does not require a detailed model of 
the heating grid. It can be parameterized from measurement data and hence is ideal 
for brownfield installations. Only load share of and transport time to consumption 
areas are needed. The delay matrix approach enables very fast solution times as it 
assumes constant transport delays except for the extension to multiple CHPs. Here, 
additional binaries representing the different velocities are introduced, enabling a 
variable delay at the cost of slightly increased solution times. Nevertheless, compared 
to the other algorithms developed in this thesis it is still very fast as shown in the case 
study in Section 6.1. In Section 6.2 the delay matrix approach was applied to a real 
world heating grid. To reduce modeling efforts the delay matrix was parameterized 
using measurement data. The results show that the delay matrix approach can reduce 
operational cost by up to 6 % or several thousand Euros per day which could easily 
result in in savings of more than 180,000 € for one winter season. The variations of the 
savings are dependent on optimization horizon (the longer the better), grid topology 
(a few larger pipes are better than many small pipes) and ambient air temperature 
(at very low temperatures the heating grid is fully loaded and thermal storage is 
not possible). In addition, the solutions in the real world case study were almost all 
achieved in less than 1 s. Only few optimization runs took up to 100 s. Thus, the 
delay matrix approach is suitable for every day operations planning. 


Taken together, this thesis hopefully enables a wider usage of scheduling of CHPs 
considering heating grid dynamics in heating grids worldwide to leverage its full 
potential of economic and environmental benefits. Besides direct application of the 
developed approaches, it might help researchers to improve their algorithms with a 
globally optimal benchmark. 
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